OFFSET
1,1
COMMENTS
Indices of primes in A066606.
LINKS
Michael S. Branicky, Table of n, a(n) for n = 1..10000 (terms 1..1000 from Alois P. Heinz)
MAPLE
a:= proc(n) option remember; local k;
for k from 1+ `if`(n=1, 0, a(n-1))
while not isprime(2 &^k mod k^2) do od; k
end:
seq (a(n), n=1..100); # Alois P. Heinz, Nov 17 2012
MATHEMATICA
Flatten[Position[Table[PowerMod[2, k, k^2], {k, 1000}], _?(PrimeQ[#] &)]] (* T. D. Noe, Nov 15 2012 *)
Select[Range[1000], PrimeQ[PowerMod[2, #, #^2]]&] (* Harvey P. Dale, Mar 29 2020 *)
PROG
(Java)
import java.math.BigInteger;
public class A219219 {
public static void main (String[] args) {
BigInteger b2 = BigInteger.valueOf(2);
for (int n=1; ; n++) {
BigInteger bn = BigInteger.valueOf(n);
BigInteger pp = b2.modPow(bn, bn.multiply(bn));
if (pp.isProbablePrime(2)) {
if (pp.isProbablePrime(80))
System.out.printf("%d, ", n);
}
}
}
}
(Python)
from sympy import isprime
def aupto(limit):
alst = []
for k in range(1, limit+1):
if isprime(pow(2, k, k*k)): alst.append(k)
return alst
print(aupto(939)) # Michael S. Branicky, May 21 2021
CROSSREFS
KEYWORD
nonn
AUTHOR
Alex Ratushnyak, Nov 15 2012
STATUS
approved