%I #38 Jul 28 2021 04:26:22
%S 5,13,563,1277,780887
%N Numbers k such that (k-1)^2 < (k-1)! mod k^2.
%C The Wilson primes (A007540) are terms of this sequence.
%C a(n) is prime or twice a prime. Otherwise (k-1)! mod k^2 = 0 for k > 9 where k is not a prime and not twice a prime. - _David A. Corneth_, Jul 23 2017
%t Select[Range[10^4],(#-1)^2<Mod[(#-1)!,#^2]&] (* _Giorgos Kalogeropoulos_, Jul 23 2021 *)
%o (PARI) for(n=1,1e5,a=(n-1)!%n^2;if((n-1)^2<a,print1(n", ")))
%o (PARI) is(n) = (n-1)^2 < lift(Mod((n-1)!, n^2)) \\ _Felix Fröhlich_, Jul 23 2017
%o (PARI) val(n, p) = my(r=0); while(n, r+=n\=p); r
%o is(n) = my(f = factor(n), r = Mod(1, n^2)); if(#f~ > 2, return(0), if(#f~==2, if(f[1,1]!=2, return(0)))); forprime(p=2,n-1, r*=Mod(p,n^2)^val(n-1,p)); (n-1)^2 < lift(r) \\ _David A. Corneth_, Jul 23 2017
%o (Python)
%o def ok(n):
%o nn = n**2; f = 1%nn
%o for k in range(1, n): f = f*k%nn
%o return (n-1)**2 < f
%o print(list(filter(ok, range(1, 1300)))) # _Michael S. Branicky_, Jul 23 2021
%o (Python) # faster for initial segment of sequence
%o from math import factorial
%o def afind(limit, startk=1):
%o k = startk; kkprev = (k-1)**2; f = factorial(k-1)
%o while k < limit:
%o kk = k*k
%o if kkprev < f%kk: print(k, end=", ")
%o kkprev = kk; f *= k; k += 1
%o afind(10000) # _Michael S. Branicky_, Jul 25 2021
%Y Cf. A007540.
%K nonn,hard,more
%O 1,1
%A _Gionata Neri_, Jul 23 2017
%E a(5) from _Chai Wah Wu_, Jul 30 2017