%I #22 Nov 17 2019 09:33:19
%S 9,91,1729,4187,6533,8149,8401,10001,11111,19201,21931,50851,79003,
%T 83119,94139,100001,102173,118301,118957,134863,139231,148417,158497,
%U 166499,188191,196651,201917,216001,226273,231337,237169,251251,287809,302177
%N Strong pseudoprimes to base 10.
%H Amiram Eldar, <a href="/A020236/b020236.txt">Table of n, a(n) for n = 1..10000</a> (terms 1..203 from R. J. Mathar)
%H <a href="/index/Ps#pseudoprimes">Index entries for sequences related to pseudoprimes</a>
%e From _Alonso del Arte_, Aug 10 2018: (Start)
%e 9 is a strong pseudoprime to base 10. It's not enough to check that 10^8 = 1 mod 9. Since 8 = 1 * 2^3, we also need to verify that 10 = 1 mod 9 and 10^2 = 1 mod 9 as well. Since these are both equal to 1, we see that 9 is indeed a strong pseudoprime to base 10.
%e 91 is also a strong pseudoprime to base 10. Besides checking that 10^90 = 1 mod 91, since 90 = 45 * 2, we also check that 10^45 = -1 mod 91; the -1 is enough to satisfy the definition of a strong pseudoprime.
%e 99 is a Fermat pseudoprime to base 10 (see A005939) but it is not a strong pseudoprime to base 10. Although 10^98 = 1 mod 99, since 98 = 49 * 2, we have to check 10^49 mod 99, and there we find not -1 nor 1 but 10. Therefore 99 is not in this sequence. (End)
%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[10, #] &] (* _Alonso del Arte_, Aug 10 2018 *)
%Y Cf. A005939, A001262, A020229, A020230, A020231, A020232, A020233, A020234, A020235.
%K nonn
%O 1,1
%A _David W. Wilson_