A219219 Numbers k such that 2^k (mod k^2) is prime.

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

Author

Alex Ratushnyak, Nov 15 2012

Keywords

• nonn

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)

Crossrefs

Cf. A066606.

Programs

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);
      }
    }
  }
}

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 *)

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