Reminder: The OEIS is hiring a new managing editor, and the application deadline is January 26.
%I #18 May 08 2024 02:28:00
%S 3,3,53,503,4297,947,10589,17903,624401,7151083,45543077,30611047,
%T 612126937,2280521251,649288301,26566080479,28921314337,303937208923,
%U 1086758949557,12299159511127,39118361784041,18314722943123,64249761922429,2484777068103119,1148475719438129,14810825716436683
%N a(n) is the least prime p such that p^(2*n+1) == 2*n+1 (mod 2^(2*n+1)).
%H Robert Israel, <a href="/A339758/b339758.txt">Table of n, a(n) for n = 0..1654</a>
%e For n = 2, 2*n+1 = 5, and 53 is the least prime q such that q^5 == 5 (mod 2^5), so a(2) = 53.
%p f:= proc(k) local x,m;
%p for m from subs(msolve(x^k=k, 2^k),x) by 2^k do
%p if isprime(m) then return m fi
%p od
%p end proc:
%p seq(f(2*i+1),i=0..50);
%o (PARI) a(n) = my(p=2); while (Mod(p, 2^(2*n+1))^(2*n+1) != 2*n+1, p = nextprime(p+1)); p; \\ _Michel Marcus_, Dec 16 2020
%o (Python)
%o from itertools import count
%o from sympy import nthroot_mod, isprime
%o def A339758(n):
%o m = (n<<1)+1
%o r = 1<<m
%o a = sorted(nthroot_mod(m,m,r,all_roots=True))
%o for i in count(0):
%o for k in a:
%o if isprime(k+i*r):
%o return int(k+i*r) # _Chai Wah Wu_, May 07 2024
%K nonn
%O 0,1
%A _J. M. Bergot_ and _Robert Israel_, Dec 16 2020