login
Strong pseudoprimes to base 54.
1

%I #30 Aug 06 2018 07:26:22

%S 55,341,361,385,1247,2701,2863,4033,4069,7957,9073,14905,15409,21349,

%T 22495,27157,29341,32689,33227,37921,42001,42121,42127,49141,55831,

%U 56449,60701,62893,70801,77293,83333,107929,128143,132193,145921,150553,152497

%N Strong pseudoprimes to base 54.

%H R. J. Mathar, <a href="/A020280/b020280.txt">Table of n, a(n) for n = 1..300</a>

%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/StrongPseudoprime.html">Strong Pseudoprime</a>

%H <a href="/index/Ps#pseudoprimes">Index entries for sequences related to pseudoprimes</a>

%t strongPseudoprimeQ[b_, n_] := Module[{rems = Table[PowerMod[b, (n - 1)/2^expo, n], {expo, 0, IntegerExponent[n - 1, 2]}]}, (rems[[-1]] == 1 || MemberQ[rems, n - 1]) && PowerMod[b, n - 1, n] == 1]; max = 5000; Select[Complement[Range[2, max], Prime[Range[PrimePi[max]]]], strongPseudoprimeQ[54, #] &] (* _Alonso del Arte_, Aug 03 2018 *)

%o (PARI)

%o oddpart(n) = if (n % 2, n, oddpart(n/2))

%o isA020280(n) = {local(d, s, res); d = oddpart(n - 1); s = bigomega((n - 1)/d); /* count factors of 2 */

%o res = 0; if(s != 0 & !isprime(n), /* n is odd and composite */

%o if (Mod(54, n)^d == Mod(1, n), res = 1, /* 54^d = 1 (mod n) */

%o for (r = 0, s - 1, if (Mod(54, n)^(d * 2^r) == Mod(-1, n), res = 1)))); res} // 54^(d * 2^r) = -1 (mod n) \\ _Michael B. Porter_, Oct 23 2009

%K nonn

%O 1,1

%A _David W. Wilson_