#!/usr/bin/env python3
from fractions import Fraction
import json, math, sys
from pathlib import Path

ZERO=(Fraction(0),Fraction(0))
ONE=(Fraction(1),Fraction(0))
PHI=(Fraction(1,2),Fraction(1,2))

def F(s): return Fraction(s)
def pair(o): return (F(o["a"]),F(o["b"]))
def add(u,v): return (u[0]+v[0],u[1]+v[1])
def sub(u,v): return (u[0]-v[0],u[1]-v[1])
def mul(u,v): return (u[0]*v[0]+5*u[1]*v[1],u[0]*v[1]+u[1]*v[0])
def inv(u):
    d=u[0]*u[0]-5*u[1]*u[1]
    assert d
    return (u[0]/d,-u[1]/d)
def div(u,v): return mul(u,inv(v))
def power(u,n):
    a=ONE
    while n:
        if n&1: a=mul(a,u)
        u=mul(u,u); n//=2
    return a
def scale(u,n): return (u[0]*n,u[1]*n)
def sign(u):
    a,b=u
    if b==0: return (a>0)-(a<0)
    if a==0: return (b>0)-(b<0)
    if a>0 and b>0: return 1
    if a<0 and b<0: return -1
    d=a*a-5*b*b
    if d==0: return 0
    if a>0: return 1 if d>0 else -1
    return 1 if d<0 else -1

def pmul(A,B):
    C={}
    for ma,ca in A.items():
        for mb,cb in B.items():
            m=(ma[0]+mb[0],ma[1]+mb[1],ma[2]+mb[2])
            C[m]=add(C.get(m,ZERO),mul(ca,cb))
    return {m:c for m,c in C.items() if c!=ZERO}

def ppow(P,n):
    R={(0,0,0):ONE}
    while n:
        if n&1: R=pmul(R,P)
        n//=2
        if n: P=pmul(P,P)
    return R

def dfodd(n):
    if n<=0: return 1
    r=1
    for k in range(n,0,-2): r*=k
    return r

def sphere_mean(P):
    out=ZERO
    for (i,j,k),c in P.items():
        if i%2 or j%2 or k%2: continue
        q=Fraction(dfodd(i-1)*dfodd(j-1)*dfodd(k-1),dfodd(i+j+k+1))
        out=add(out,scale(c,q))
    return out

def eval_poly(P,xyz):
    caches=[]
    maxdeg=max(max(m) for m in P)
    for q in xyz:
        c=[ONE]
        for _ in range(maxdeg): c.append(mul(c[-1],q))
        caches.append(c)
    out=ZERO
    for (i,j,k),coef in P.items():
        out=add(out,mul(coef,mul(caches[0][i],mul(caches[1][j],caches[2][k]))))
    return out

def inv_values(part,H6,H10):
    a,b,c=part
    px=(Fraction(c-b),Fraction(0))
    py=add((Fraction(a),Fraction(0)),scale(PHI,b+c))
    pz=scale(PHI,a)
    n2=add(mul(px,px),add(mul(py,py),mul(pz,pz)))
    u6=div(eval_poly(H6,(px,py,pz)),power(n2,3))
    u10=div(eval_poly(H10,(px,py,pz)),power(n2,5))
    return u6,u10

def monoms(t):
    q=[]
    for r in range(t//6+1):
        for s in range(t//10+1):
            d=6*r+10*s
            if d<=t: q.append((d,r,s))
    return [(r,s) for d,r,s in sorted(q)]

def main():
    path=Path(sys.argv[1] if len(sys.argv)>1 else "classI_exact_certificates_f2_f21.json")
    D=json.loads(path.read_text())
    H6={tuple(e["monomial"]):pair(e["coeff"]) for e in D["H6"]}
    H10={tuple(e["monomial"]):pair(e["coeff"]) for e in D["H10"]}
    max_tau=max(v["tau"] for v in D["frequencies"].values())
    needed=monoms(max_tau)
    P6={0:{(0,0,0):ONE}}
    P10={0:{(0,0,0):ONE}}
    for r in range(1,max(r for r,s in needed)+1): P6[r]=pmul(P6[r-1],H6)
    for s in range(1,max(s for r,s in needed)+1): P10[s]=pmul(P10[s-1],H10)
    rhs={(r,s):sphere_mean(pmul(P6[r],P10[s])) for r,s in needed}

    total_eq=0
    for fs,rec in sorted(D["frequencies"].items(),key=lambda kv:int(kv[0])):
        f=int(fs); tau=int(rec["tau"])
        parts=[tuple(w["partition"]) for w in rec["weights"]]
        W=[pair(w["weight"]) for w in rec["weights"]]
        assert all(sign(w)>0 for w in W), f"nonpositive weight at f={f}"
        assert sum((w[0] for w in W),Fraction(0))==1 and sum((w[1] for w in W),Fraction(0))==0, f"sum fail f={f}"
        UV=[inv_values(p,H6,H10) for p in parts]
        for r,s in monoms(tau):
            lhs=ZERO
            for w,(u6,u10) in zip(W,UV):
                lhs=add(lhs,mul(w,mul(power(u6,r),power(u10,s))))
            assert lhs==rhs[(r,s)], f"moment fail f={f}, (r,s)=({r},{s})"
            total_eq+=1
        print(f"PASS f={f:2d} tau={tau:2d} orbits={len(W):2d} minW={min(float(w[0])+float(w[1])*math.sqrt(5) for w in W):.12g}")
    print(f"ALL PASS: {len(D['frequencies'])} frequencies, {total_eq} invariant moment equations, exact Q(sqrt(5)) arithmetic.")

if __name__=="__main__":
    main()