"""Machine checks for the full cyclic-rigidity Newton-edge proof.

The proof itself is symbolic/algebraic; this file verifies:
  1. the exact cyclic Keller determinant identity;
  2. the general Newton-edge differential identity with coprime p,q;
  3. the nonzero terminal coefficient obstruction;
  4. bounded exact linear-algebra searches for arbitrary multi-mode g.
No floating point is used.
"""
from __future__ import annotations
import hashlib, json, platform, random
from pathlib import Path
import sympy as sp

u,y,x,t=sp.symbols('u y x t')
m,p,q,r,s,M,K=sp.symbols('m p q r s M K', positive=True, integer=True)
Phi=sp.Function('Phi')(t); A=sp.Function('A')(t)

# 1: direct determinant identity on exact random polynomials.
def random_poly(maxu,maxy,rng):
    out=0
    for i in range(maxu+1):
      for j in range(maxy+1):
        z=rng.randint(-3,3)
        if z: out += sp.Rational(z,rng.randint(1,4))*u**i*y**j
    return sp.expand(out)

rng=random.Random(20260721)
identity_checks=0
for mi in range(1,8):
  for _ in range(4):
    H=random_poly(3,4,rng); g=random_poly(3,3,rng)
    F1=H.subs(u,x**mi); F2=x*g.subs(u,x**mi)
    direct=sp.expand(sp.det(sp.Matrix([[sp.diff(F1,x),sp.diff(F1,y)],
                                      [sp.diff(F2,x),sp.diff(F2,y)]])))
    reduced=sp.expand((mi*u*sp.diff(H,u)*sp.diff(g,y)
                       -sp.diff(H,y)*(g+mi*u*sp.diff(g,u))).subs(u,x**mi))
    assert sp.expand(direct-reduced)==0
    identity_checks+=1

# 2: edge identity, treating t=u^q*y^p by explicit chain rule.
# H_edge=u^r*y^s*Phi(t), g_edge=A(t).
# D_edge/u^r/y^(s-1) must equal the displayed bracket.
Hu_factor = r*Phi + q*t*sp.diff(Phi,t)         # u*H_u /(u^r y^s)
Hy_factor = s*Phi + p*t*sp.diff(Phi,t)         # H_y /(u^r y^(s-1))
gy_factor = p*t*sp.diff(A,t)                   # y*g_y
ugu_factor = q*t*sp.diff(A,t)                  # u*g_u
edge_got=sp.expand(m*Hu_factor*gy_factor-Hy_factor*(A+m*ugu_factor))
sigma=(q*s-p*r)/q
edge_want=sp.expand(-p*t*A*sp.diff(Phi,t)-s*A*Phi-m*q*sigma*t*sp.diff(A,t)*Phi)
assert sp.simplify(edge_got-edge_want)==0

# 3: terminal coefficient is strictly nonzero for actual positive integers.
terminal_checks=0
for mi in range(1,10):
 for pi in range(1,8):
  for qi in range(1,8):
   if sp.gcd(pi,qi)!=1: continue
   for ri in range(0,qi):
    for si in range(1,6):
     sig=sp.Rational(qi*si-pi*ri,qi)
     if sig<=0: continue
     for Mi in range(0,5):
      for Ki in range(1,5):
       coeff=pi*Mi+si+mi*qi*sig*Ki
       assert coeff>0
       terminal_checks+=1

# 4: bounded modular exclusions. Inconsistency modulo a prime certifies inconsistency over Q.
PRIME=1000003
def mons(U,Y): return [(i,j) for i in range(U+1) for j in range(Y+1)]
def dict_poly(expr): return sp.Poly(sp.expand(expr),u,y).as_dict()
def no_slice_mod_prime(mi,g,U,Y,prime=PRIME):
    columns=[]
    for i,j in mons(U,Y):
      h=u**i*y**j
      z=dict_poly(mi*u*sp.diff(h,u)*sp.diff(g,y)-sp.diff(h,y)*(g+mi*u*sp.diff(g,u)))
      columns.append({k:int(v)%prime for k,v in z.items() if int(v)%prime})
    basis={}
    for v in columns:
      v=dict(v)
      while v:
       pivot=max(v)
       if pivot not in basis:
        inv=pow(v[pivot],prime-2,prime)
        basis[pivot]={k:(z*inv)%prime for k,z in v.items() if z%prime}; break
       f=v[pivot]
       for k,z in basis[pivot].items():
        v[k]=(v.get(k,0)-f*z)%prime
        if not v[k]: v.pop(k,None)
    target={(0,0):1}
    while target:
      pivot=max(target)
      if pivot not in basis:return True
      f=target[pivot]
      for k,z in basis[pivot].items():
       target[k]=(target.get(k,0)-f*z)%prime
       if not target[k]:target.pop(k,None)
    return False

bounded=[]
for mi in range(1,7):
 for case in range(20):
    g=sp.Integer(rng.randint(1,5))
    for i in range(1,rng.randint(2,4)+1):
      ai=sum(sp.Integer(rng.randint(-3,3))*y**j for j in range(rng.randint(0,4)+1))
      g += u**i*ai
    g=sp.expand(g)
    if sp.degree(g,u)==0: continue
    excluded=no_slice_mod_prime(mi,g,7,9)
    assert excluded,(mi,g)
    bounded.append({'m':mi,'g':str(g),'H_u_cap':7,'H_y_cap':9,'excluded_mod_prime':PRIME})

cert={
 'status':'PASS',
 'theorem':'full cyclic rigidity Newton-edge checks',
 'exact_arithmetic':True,
 'determinant_identity_checks':identity_checks,
 'edge_identity':'PASS',
 'terminal_coefficient_positive_checks':terminal_checks,
 'bounded_multimode_exclusions':len(bounded),
 'bounded_caps':{'u':7,'y':9,'prime':PRIME},
 'sympy':sp.__version__,
 'python':platform.python_version(),
 'edge_formula':'-p*t*A*Phi_prime - s*A*Phi - m*q*sigma*t*A_prime*Phi',
 'terminal_coefficient':'p*deg(Phi)+s+m*q*sigma*deg(A) > 0'
}
out=Path(__file__).with_name('full_cyclic_rigidity_certificate.json')
out.write_text(json.dumps(cert,indent=2)+'\n',encoding='utf-8')
print(json.dumps(cert,indent=2))
print('certificate_sha256',hashlib.sha256(out.read_bytes()).hexdigest())