%I #17 Aug 03 2023 11:06:54
%S 5,13,17,29,41,53,73,97,137,229,241,257,281,313,397,409,449,457,461,
%T 521,541,617,653,661,677,733,769,809,853,857,929,953,997,1021,1069,
%U 1109,1201,1213,1217,1249,1277,1361,1373,1409,1489,1553,1597,1609,1621,1697
%N Lagrange primes: primes p == 1 (mod 4) such that X = (p-1)/2 is the least solution in the interval [1, (p-1)/2] of the congruence (X!)^2 == -1 (mod p).
%D J. B. Cosgrave, A Mersenne-Wieferich Odyssey, Manuscript, May 2022. See Section 18.5.
%H Michael S. Branicky, <a href="/A354155/b354155.txt">Table of n, a(n) for n = 1..3100</a>
%o (Python)
%o from itertools import islice
%o from sympy import factorial, nextprime
%o def agen(): # generator of terms
%o p = 5
%o while True:
%o X = (p-1)//2
%o Xf = factorial(X)**2
%o if all(pow(factorial(Y), 2, p)+1 != p for Y in range(X-1, 0, -1)):
%o yield p
%o p = nextprime(p)
%o while p%4 != 1:
%o p = nextprime(p)
%o print(list(islice(agen(), 5))) # _Michael S. Branicky_, May 30 2022
%o (PARI) is(n)=if(n%4 != 1 || !isprime(n), return(0)); my(t1=lift(sqrt(Mod(-1,n))), t2=n-t1, t=Mod(1,n)); for(k=2,n\2, if(t==t1 || t==t2, return(0)); t*=k); 1 \\ _Charles R Greathouse IV_, Aug 03 2023
%o (PARI) list(lim)=my(v=List()); forprimestep(p=5,lim\1,4, my(t1=lift(sqrt(Mod(-1,p))),t2=p-t1, t=Mod(1,p)); for(k=2,p\2, if(t==t1 || t==t2, next(2)); t*=k); listput(v,p)); Vec(v) \\ _Charles R Greathouse IV_, Aug 03 2023
%Y Cf. A002144, A354156.
%K nonn
%O 1,1
%A _N. J. A. Sloane_, May 29 2022, based on Section 18.5 of Cosgrave (2022)
%E a(20) and beyond from _Michael S. Branicky_, May 30 2022