login

Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).

A219219
Numbers k such that 2^k (mod k^2) is prime.
1
5, 21, 53, 55, 61, 95, 111, 155, 165, 189, 193, 213, 221, 227, 245, 249, 257, 289, 291, 303, 305, 307, 317, 339, 345, 355, 363, 383, 385, 423, 429, 437, 457, 465, 477, 505, 577, 597, 601, 607, 621, 653, 655, 679, 705, 715, 727, 749, 751, 765, 781, 849, 889, 939
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
Cf. A066606.
Sequence in context: A316435 A339623 A272013 * A272810 A147834 A160378
KEYWORD
nonn
AUTHOR
Alex Ratushnyak, Nov 15 2012
STATUS
approved