|
%I #30 Nov 21 2019 04:18:31
%S 2,4,6,2,10,12,2,16,18,2,22,4,2,28,30,2,2,36,2,40,42,4,46,6,2,52,2,2,
%T 58,60,2,8,66,2,70,72,2,2,78,2,82,8,2,88,18,2,2,96,2,100,102,8,106,
%U 108,2,112,2,4,2,10,2,4,126,2,130,18,2,136,138,2,2,8,2
%N Number of witnesses for Solovay-Strassen primality test of 2*n+1.
%C Number of bases b, 1 <= b <= 2*n, such that GCD(b, 2*n+1) = 1 and b^n == (b / 2*n+1) (mod 2*n+1), where (b / 2*n+1) is a Jacobi symbol.
%C If 2*n+1 is composite then it is the number of bases b, 1 <= b <= 2*n, in which 2*n+1 is an Euler-Jacobi pseudoprime.
%C Differs from A071294 from n = 22.
%D Paulo Ribenboim, The Little Book of Bigger Primes, 2nd ed., Springer-Verlag, New York, 2004, p. 96.
%H Amiram Eldar, <a href="/A329726/b329726.txt">Table of n, a(n) for n = 1..10000</a>
%H Louis Monier, <a href="https://doi.org/10.1016/0304-3975(80)90007-9">Evaluation and comparison of two efficient primality testing algorithms</a>, Theoretical Computer Science, Vol. 11 (1980), pp. 97-108.
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/Euler-JacobiPseudoprime.html">Euler-Jacobi Pseudoprime</a>.
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Euler-Jacobi_pseudoprime">Euler-Jacobi pseudoprime</a>.
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Solovay%E2%80%93Strassen_primality_test">Solovay-Strassen primality test</a>.
%F a(n) = delta(n) * Product_{p|n} gcd((n-1)/2, p-1), where delta(n) = 2 if nu(n-1, 2) = min_{p|n} nu(p-1, 2), 1/2 if there is a prime p|n such that nu(p, n) is odd and nu(p-1, 2) < nu(n-1, 2), and 1 otherwise, where nu(n, p) is the exponent of the highest power of p dividing n.
%F a(p) = p-1 for prime p.
%e a(1) = 2 since there are 2 bases b in which 2*1 + 1 = 3 is an Euler-Jacobi pseudoprime: b = 1 since GCD(1, 3) = 1 and 1^1 == (1 / 3) == 1 (mod 3), and b = 2 since GCD(2, 3) = 1 and 2^1 == (2 / 3) == -1 (mod 3).
%t v[n_] := Min[IntegerExponent[#, 2]& /@ (FactorInteger[n][[;;, 1]] - 1)];
%t pQ[n_, p_] := OddQ[IntegerExponent[n, p]] && IntegerExponent[p-1, 2] < IntegerExponent[n-1, 2];
%t psQ[n_] := AnyTrue[FactorInteger[n][[;;, 1]], pQ[n, #] &];
%t delta[n_] := If[IntegerExponent[n-1, 2] == v[n], 2, If[psQ[n], 1/2, 1]];
%t a[n_] := delta[n] * Module[{p = FactorInteger[n][[;;, 1]]}, Product[GCD[(n-1)/2, p[[k]]-1], {k, 1, Length[p]}]];
%t Table[a[n], {n, 3, 147, 2}]
%Y Cf. A006093, A007324, A007814, A047713, A063994, A071294.
%K nonn
%O 1,1
%A _Amiram Eldar_, Nov 20 2019
|
