#!/usr/bin/env python3
import json,sys,hashlib,math
from decimal import Decimal, getcontext
from fractions import Fraction as F
from math import factorial,sqrt
sys.set_int_max_str_digits(1000000)

class K5:
 __slots__=('a','b')
 def __init__(self,a=0,b=0):self.a=a if isinstance(a,F) else F(a);self.b=b if isinstance(b,F) else F(b)
 def __add__(self,o):
  if not isinstance(o,K5):o=K5(o)
  return K5(self.a+o.a,self.b+o.b)
 __radd__=__add__
 def __neg__(self):return K5(-self.a,-self.b)
 def __sub__(self,o):return self+(-o)
 def __rsub__(self,o):return K5(o)-self
 def __mul__(self,o):
  if not isinstance(o,K5):o=K5(o)
  return K5(self.a*o.a+5*self.b*o.b,self.a*o.b+self.b*o.a)
 __rmul__=__mul__
 def inv(self):
  d=self.a*self.a-5*self.b*self.b
  if d==0:raise ZeroDivisionError
  return K5(self.a/d,-self.b/d)
 def __truediv__(self,o):
  if not isinstance(o,K5):o=K5(o)
  return self*o.inv()
 def __pow__(self,n):
  if n<0:return self.inv()**(-n)
  r=K5(1);b=self
  while n:
   if n&1:r=r*b
   b=b*b;n//=2
  return r
 def __eq__(self,o):
  if not isinstance(o,K5):o=K5(o)
  return self.a==o.a and self.b==o.b
 def zero(self):return self.a==0 and self.b==0
 def sign(self):
  if self.b==0:return (self.a>0)-(self.a<0)
  if self.a>=0 and self.b>=0:return 1
  if self.a<=0 and self.b<=0:return -1
  aa=self.a*self.a;bb=5*self.b*self.b
  if aa==bb:return 0
  if self.a>0:return 1 if aa>bb else -1
  return -1 if aa>bb else 1
 def decimal(self,prec=90):
  getcontext().prec=prec
  return (Decimal(self.a.numerator)/Decimal(self.a.denominator) +
          (Decimal(self.b.numerator)/Decimal(self.b.denominator))*Decimal(5).sqrt())
 def __lt__(self,o):
  if not isinstance(o,K5):o=K5(o)
  return (self-o).sign()<0

def Q(s):return F(s)
SQ=K5(0,1);PHI=(K5(1)+SQ)/K5(2);R2=PHI+K5(2)
VERT=[]
for s1 in (1,-1):
 for s2 in (1,-1):VERT += [(K5(),K5(s1),K5(s2)*PHI),(K5(s1),K5(s2)*PHI,K5()),(K5(s1)*PHI,K5(),K5(s2))]
FACE=((K5(),K5(1),PHI),(K5(-1),PHI,K5()),(K5(1),PHI,K5()))
def dot(x,y):return sum((x[i]*y[i] for i in range(3)),K5())
def padd(A,B):
 d=A.copy()
 for e,c in B.items():
  z=d.get(e,K5())+c
  if z.zero():d.pop(e,None)
  else:d[e]=z
 return d
def pscale(A,s):return {e:c*s for e,c in A.items() if not (c*s).zero()}
def pmul(A,B):
 d={}
 for (i,j,k),a in A.items():
  for (u,v,w),b in B.items():
   e=(i+u,j+v,k+w);d[e]=d.get(e,K5())+a*b
 return {e:c for e,c in d.items() if not c.zero()}
def linpow(v,n):
 d={};vx,vy,vz=v
 for i in range(n+1):
  for j in range(n-i+1):
   k=n-i-j;c=K5(F(factorial(n),factorial(i)*factorial(j)*factorial(k)))*vx**i*vy**j*vz**k
   if not c.zero():d[(i,j,k)]=c
 return d
def invpoly(n):
 d={}
 for v in VERT:d=padd(d,linpow(v,n))
 return pscale(d,R2**(-(n//2)))
U6=invpoly(6);U10=invpoly(10)
def odddf(n):
 r=1
 for x in range(n,0,-2):r*=x
 return r
def sphmono(i,j,k):
 if i%2 or j%2 or k%2:return F(0)
 a,b,c=i//2,j//2,k//2
 return F(odddf(2*a-1)*odddf(2*b-1)*odddf(2*c-1),odddf(2*(a+b+c)+1))
def sphint(P):
 r=K5()
 for (i,j,k),c in P.items():r+=c*K5(sphmono(i,j,k))
 return r
P6={0:{(0,0,0):K5(1)},1:U6};P10={0:{(0,0,0):K5(1)},1:U10}
def moment(r,s):
 for n in range(2,r+1):
  if n not in P6:P6[n]=pmul(P6[n-1],U6)
 for n in range(2,s+1):
  if n not in P10:P10[n]=pmul(P10[n-1],U10)
 return sphint(pmul(P6[r],P10[s]))
def pairs(t):
 z=[]
 for s in range(t//10+1):
  for r in range((t-10*s)//6+1):z.append((r,s))
 return sorted(z,key=lambda q:(6*q[0]+10*q[1],q[1],q[0]))
def parts(f):
 z=[]
 for a in range(f,-1,-1):
  for b in range(min(a,f-a),-1,-1):
   c=f-a-b
   if 0<=c<=b:z.append((a,b,c))
 return z
def orbit(rep):
 a,b,c=rep;p=tuple(K5(a)*FACE[0][i]+K5(b)*FACE[1][i]+K5(c)*FACE[2][i] for i in range(3));ns=dot(p,p);s6=K5();s10=K5()
 for v in VERT:
  q=dot(p,v);s6+=q**6;s10+=q**10
 return s6/(R2**3*ns**3),s10/(R2**5*ns**5)
def fracmod(q,p):return (q.numerator%p)*pow(q.denominator%p,-1,p)%p
def kmod(x,p,r):return (fracmod(x.a,p)+r*fracmod(x.b,p))%p
def detmod(M,p):
 n=len(M);d=1;sg=1
 for k in range(n):
  q=next((i for i in range(k,n) if M[i][k]),None)
  if q is None:return 0
  if q!=k:M[k],M[q]=M[q],M[k];sg=-sg
  pv=M[k][k];d=d*pv%p;inv=pow(pv,-1,p)
  for i in range(k+1,n):
   fac=M[i][k]*inv%p
   if fac:
    for j in range(k+1,n):M[i][j]=(M[i][j]-fac*M[k][j])%p
    M[i][k]=0
 return d if sg>0 else (-d)%p

def main(path):
 raw=open(path,'rb').read();D=json.loads(raw)
 assert D['format']=='CLASS_I_CUBATURE_F23_F24_EXACT/1.1'
 p=int(D['rank_witness']['prime']);rp=int(D['rank_witness']['sqrt5_root_mod_prime']);rm=int(D['rank_witness']['sqrt5_negative_root_mod_prime'])
 assert rp*rp%p==5%p and rm*rm%p==5%p and (rp+rm)%p==0
 total=0
 for fs in ('23','24'):
  f=int(fs);rec=D['frequencies'][fs];t=int(rec['target_degree']);reps=parts(f);ps=pairs(t)
  assert rec['nodes']==10*f*f+2
  assert len(reps)==rec['orbits'] and len(ps)==rec['moment_conditions']
  assert rec['rank']==len(ps) and rec['affine_dimension']==len(reps)-len(ps)
  weights=rec['weights'];assert [tuple(x['partition']) for x in weights]==reps
  W=[K5(Q(x['weight']['a']),Q(x['weight']['b'])) for x in weights]
  assert all(w.sign()>0 for w in W)
  assert sum(W,K5())==K5(1)
  vals=[orbit(rep) for rep in reps];rhs=[moment(r,s) for r,s in ps]
  A=[[(u6**r)*(u10**s) for (u6,u10) in vals] for r,s in ps]
  for ri,row in enumerate(A):
   lhs=sum((row[j]*W[j] for j in range(len(W))),K5())
   assert lhs==rhs[ri],(f,ps[ri])
   total+=1
  free={tuple(x['partition']) for x in rec['free_orbits']}
  for e in rec['free_orbits']:
   j=reps.index(tuple(e['partition']))
   assert W[j]==K5(Q(e['weight']),0)
  pivs=[j for j,r in enumerate(reps) if r not in free];assert len(pivs)==len(ps)
  dp=detmod([[kmod(A[i][j],p,rp) for j in pivs] for i in range(len(ps))],p)
  dm=detmod([[kmod(A[i][j],p,rm) for j in pivs] for i in range(len(ps))],p)
  assert dp==int(rec['rank_witness_det_plus']) and dm==int(rec['rank_witness_det_minus']) and dp and dm
  # canonical weight text hash
  lines=[f"{r[0]},{r[1]},{r[2]}:{w.a.numerator}/{w.a.denominator}:{w.b.numerator}/{w.b.denominator}" for r,w in zip(reps,W)]
  dg=hashlib.sha256(('\n'.join(lines)+'\n').encode()).hexdigest();assert dg==rec['canonical_weight_text_sha256']
  mi=0
  for i in range(1,len(W)):
   if W[i] < W[mi]: mi=i
  assert tuple(rec['minimum_orbit'])==reps[mi]
  assert rec['minimum_weight']=={'a':str(W[mi].a),'b':str(W[mi].b)}
  d60=format(W[mi].decimal(),'.60f')
  assert d60==rec['minimum_decimal_60']
  print(f"PASS f={f} tau={t} nodes={10*f*f+2} orbits={len(reps)} rows={len(ps)} rank={len(ps)} free={len(reps)-len(ps)} min={reps[mi]} ~ {d60}")
 print(f"ALL PASS: {total} exact invariant moment identities; strict positivity and rank certified over Q(sqrt(5)).")
 print('certificate_sha256',hashlib.sha256(raw).hexdigest())
if __name__=='__main__':main(sys.argv[1] if len(sys.argv)>1 else 'ClassI_Cubature_f23_f24_exact_certificate.json')