Reminder: The OEIS is hiring a new managing editor, and the application deadline is January 26.
%I #17 May 21 2021 11:28:48
%S 5,21,53,55,61,95,111,155,165,189,193,213,221,227,245,249,257,289,291,
%T 303,305,307,317,339,345,355,363,383,385,423,429,437,457,465,477,505,
%U 577,597,601,607,621,653,655,679,705,715,727,749,751,765,781,849,889,939
%N Numbers k such that 2^k (mod k^2) is prime.
%C Indices of primes in A066606.
%H Michael S. Branicky, <a href="/A219219/b219219.txt">Table of n, a(n) for n = 1..10000</a> (terms 1..1000 from Alois P. Heinz)
%p a:= proc(n) option remember; local k;
%p for k from 1+ `if`(n=1, 0, a(n-1))
%p while not isprime(2 &^k mod k^2) do od; k
%p end:
%p seq (a(n), n=1..100); # _Alois P. Heinz_, Nov 17 2012
%t Flatten[Position[Table[PowerMod[2, k, k^2], {k, 1000}], _?(PrimeQ[#] &)]] (* _T. D. Noe_, Nov 15 2012 *)
%t Select[Range[1000],PrimeQ[PowerMod[2,#,#^2]]&] (* _Harvey P. Dale_, Mar 29 2020 *)
%o (Java)
%o import java.math.BigInteger;
%o public class A219219 {
%o public static void main (String[] args) {
%o BigInteger b2 = BigInteger.valueOf(2);
%o for (int n=1; ; n++) {
%o BigInteger bn = BigInteger.valueOf(n);
%o BigInteger pp = b2.modPow(bn, bn.multiply(bn));
%o if (pp.isProbablePrime(2)) {
%o if (pp.isProbablePrime(80))
%o System.out.printf("%d, ",n);
%o }
%o }
%o }
%o }
%o (Python)
%o from sympy import isprime
%o def aupto(limit):
%o alst = []
%o for k in range(1, limit+1):
%o if isprime(pow(2, k, k*k)): alst.append(k)
%o return alst
%o print(aupto(939)) # _Michael S. Branicky_, May 21 2021
%Y Cf. A066606.
%K nonn
%O 1,1
%A _Alex Ratushnyak_, Nov 15 2012