OFFSET
1,1
COMMENTS
The number 2^k-1 is in the sequence iff k is in A054723 or in A001567. - Thomas Ordowski, Sep 02 2016
The number (2^k+1)/3 is in the sequence iff k is in A127956. - Davide Rotondo, Aug 13 2021
REFERENCES
R. K. Guy, Unsolved Problems Theory Numbers, A12.
P. Ribenboim, The Book of Prime Number Records. Springer-Verlag, NY, 2nd ed., 1989, p. 95.
LINKS
T. D. Noe, Table of n, a(n) for n = 1..10000 (using data from A001567)
Joerg Arndt, Matters Computational (The Fxtbook), section 39.10, pp. 786-792.
Chris Caldwell, Strong probable prime.
Eric Weisstein's World of Mathematics, Strong Pseudoprime.
OEIS Wiki, Strong Pseudoprime.
Wikipedia, Strong pseudoprime.
EXAMPLE
From Michael B. Porter, Sep 04 2016: (Start)
For k = 577, k-1 = 576 = 9*2^6. Since 2^(9*2^3) = 2^72 == -1 (mod 577), 577 passes the primality test, but since it is actually prime, it is not in the sequence.
For k = 3277, k-1 = 3276 = 819*2^2, and 2^(819*2) == -1 (mod 3277), so k passes the primality test, and k = 3277 = 29*113 is composite, so 3277 is in the sequence. (End)
MAPLE
A007814 := proc(n) padic[ordp](n, 2) ; end proc:
isStrongPsp := proc(n, b) local d, s, r; if type(n, 'even') or n<=1 then return false; elif isprime(n) then return false; else s := A007814(n-1) ; d := (n-1)/2^s ; if modp(b &^ d, n) = 1 then return true; else for r from 0 to s-1 do if modp(b &^ d, n) = n-1 then return true; end if; d := 2*d ; end do: return false; end if; end if; end proc:
isA001262 := proc(n) isStrongPsp(n, 2) ; end proc:
for n from 1 by 2 do if isA001262(n) then print(n); end if; end do:
# R. J. Mathar, Apr 05 2011
MATHEMATICA
sppQ[n_?EvenQ, _] := False; sppQ[n_?PrimeQ, _] := False; sppQ[n_, b_] := (s = IntegerExponent[n-1, 2]; d = (n-1)/2^s; If[PowerMod[b, d, n] == 1, Return[True], Do[If[PowerMod[b, d, n] == n-1, Return[True]]; d = 2*d, {s}]]); lst = {}; k = 3; While[k < 500000, If[sppQ[k, 2], Print[k]; AppendTo[lst, k]]; k += 2]; lst (* Jean-François Alcover, Oct 20 2011, after R. J. Mathar *)
PROG
(PARI)
isStrongPsp(n, b)={
my(s, d, r, bm) ;
if( (n% 2) ==0 || n <=1, return(0) ; ) ;
if(isprime(n), return(0) ; ) ;
s = valuation(n-1, 2) ;
d = (n-1)/2^s ;
bm = Mod(b, n)^d ;
if ( bm == Mod(1, n), return(1) ; ) ;
for(r=0, s-1,
bm = Mod(b, n)^d ;
if ( bm == Mod(-1, n),
return(1) ;
) ;
d *= 2;
) ;
return(0);
}
isA001262(n)={
isStrongPsp(n, 2)
}
{
for(n=1, 10000000000,
if(isA001262(n),
print(n)
) ;
) ;
} \\ R. J. Mathar, Mar 07 2012
(PARI) is_A001262(n, a=2)={ (bittest(n, 0) && !isprime(n) && n>8) || return; my(s=valuation(n-1, 2)); if(1==a=Mod(a, n)^(n>>s), return(1)); while(a!=-1 && s--, a=a^2); a==-1} \\ M. F. Hasler, Aug 16 2012
CROSSREFS
KEYWORD
nonn,nice,changed
AUTHOR
EXTENSIONS
More terms from David W. Wilson, Aug 15 1996
STATUS
approved