"""Reconstruct the Freudenthal divisor certificate from the displayed seed.
Generation uses SymPy and SciPy; every returned certificate is checked over Q.
No source coordinates or frozen characters are numerically specialized.
"""
import sympy as s
from scipy.optimize import linprog
from fractions import Fraction as F
from itertools import combinations
from pathlib import Path
import json,time
v=s.symbols('v1:21'); n=20

def norm_expr(expr):
    expr=s.expand(expr)
    terms=[]
    for t in expr.as_ordered_terms():
        d=t.as_powers_dict(); e=tuple(int(d.get(x,0)) for x in v)
        c=s.cancel(t/s.prod(x**a for x,a in zip(v,e)))
        terms.append((e,c))
    ell=tuple(min(e[i] for e,c in terms) for i in range(n))
    P=s.expand(expr/s.prod(x**a for x,a in zip(v,ell)))
    return ell,P

def key(expr):
    _,P=norm_expr(expr); p=s.Poly(P,*v); _,p=p.primitive()
    if p.LC()<0:p=-p
    return tuple((e,str(c)) for e,c in p.terms())

def fromkey(k):return sum(s.Rational(c)*s.prod(x**a for x,a in zip(v,e)) for e,c in k)

I1=[1,5,15,18];I2=[6,7,12,13,20];I3=[11,14];I4=[16,17,19]
eta1=(v[10]+v[2])/s.prod(v[j-1] for j in I1)
eta2=(v[1]+v[7]*v[8])/s.prod(v[j-1] for j in I2)
eta3=(v[0]*v[3]+v[18]*v[2])/(v[10]*v[13])
eta4=(v[13]*v[1]+v[19]*v[3])/(v[15]*v[16]*v[18])
t8=(v[12]*v[19]+v[16])/v[7]
t9=(v[11]*v[12]*v[18]*v[5]*v[6]+v[1])/v[8]
t10=(v[11]*v[3]+v[1]*v[4])/v[9]
first=[eta1,eta2,eta3,eta4,
(v[9]*v[15]*v[16]*v[8]+v[11]*v[12]*v[19]*v[3]*v[5])/v[1],
(v[0]*v[14]*v[16]*v[17]*v[4]+v[10]*v[13])/v[2],
(v[9]*v[14]*v[15]*v[16]*v[18]+v[10]*v[13]*v[1])/v[3],t8,t9,t10]
W=[None]*18
W[0]=v[14];W[1]=v[4];W[2]=v[17]
W[3]=(v[16]+v[13]*eta1)/v[2];W[4]=eta3
W[5]=(v[9]*v[14]+v[10]*eta4)/v[3]
W[6]=v[15]*t8/(v[12]*v[19])+v[15]*v[16]*v[8]/(v[12]*v[1]*v[19])+v[5]*t10/v[1]
W[7]=v[12];W[8]=(v[18]+v[19]*eta2)/v[8]
W[9]=v[6];W[10]=v[11];W[11]=v[5];W[12]=t8;W[13]=eta2
W[14]=(v[9]*v[6]+eta4*t9)/v[1]
W[15]=v[15]/v[2]+v[17]*t10/v[3]+v[14]*v[15]*v[16]*v[17]*v[4]*eta3/(v[2]*v[3])
W[16]=eta1;W[17]=v[0]
func=[]
for f in first+W:
 f=s.expand(f)
 if f not in func:func.append(f)
print('distinct functions',len(func),flush=True)
normalized=[norm_expr(f) for f in func]
factors={}; subsets=0; factorizations=[]
for ell,P in normalized:
 terms=s.Add.make_args(P)
 for m in range(2,len(terms)+1):
  for subset in combinations(terms,m):
   _,sub=norm_expr(sum(subset));subsets+=1
   coeff,fs=s.factor_list(sub,*v)
   assert s.expand(coeff*s.prod(z**e for z,e in fs)-sub)==0
   factorizations.append({'polynomial_terms':[[list(ex),str(cc)] for ex,cc in s.Poly(sub,*v).terms()],
      'coefficient':str(coeff), 'factors':[{'terms':[[list(ex),str(cc)] for ex,cc in s.Poly(z,*v).terms()], 'power':e} for z,e in fs]})
   for z,e in fs:
    if len(s.Add.make_args(z))>1:factors[key(z)]=fromkey(key(z))
print('subset checks',subsets,'factors',len(factors),flush=True)

def coeff_terms(P):
 return [(tuple(int(j) for j in exp),int(c)) for exp,c in s.Poly(P,*v).terms() if c]

def solve(rows,rhs,eqs):
 result=linprog([-1.]*n,A_ub=[[-float(c) for c in row] for row in rows],b_ub=[-float(c) for c in rhs], A_eq=eqs or None,b_eq=[0.]*len(eqs) if eqs else None,bounds=[(0,None)]*n,method='highs')
 assert result.success,result.message
 lam=[F(float(-z)).limit_denominator(10000) for z in result.ineqlin.marginals]
 mu=[F(float(z)).limit_denominator(10000) for z in result.eqlin.marginals]
 # c = A^t*lambda + H^t*mu + lower multipliers
 co=[sum(lam[j]*rows[j][i] for j in range(len(rows)))+sum(mu[j]*eqs[j][i] for j in range(len(eqs))) for i in range(n)]
 assert all(z<=-1 for z in co),(co,result.fun)
 upper=-sum(lam[j]*rhs[j] for j in range(len(rows)))
 assert upper in (0,1),(upper,result.fun)
 return {'upper_bound':str(upper),'inequalities':[[j,str(x)] for j,x in enumerate(lam) if x], 'equations':[[j,str(x)] for j,x in enumerate(mu) if x], 'denominator_max':max([x.denominator for x in lam+mu]+[1])}

# w=0 when all monomials are nonnegative
allrows=[]
for f in func:
 for t in s.expand(f).as_ordered_terms():
  d=t.as_powers_dict(); row=tuple(int(d.get(x,0)) for x in v)
  if row not in allrows:allrows.append(row)
zero=solve(allrows,[0]*len(allrows),[])
records=[]
U=s.Symbol('_q')
for idx,(k,qpoly) in enumerate(sorted(factors.items())):
 piv=[]
 for p,x in enumerate(v):
  pol=s.Poly(qpoly,x)
  if pol.degree()==1:
   A=pol.nth(1);B=pol.nth(0)
   if len(s.Add.make_args(A))==1 and s.gcd(A,B)==1:piv.append((p,A,B))
 assert piv, str(qpoly)
 p,A,B=piv[0];aexp=tuple(int(A.as_powers_dict().get(x,0)) for x in v)
 exps=[e for e,c in coeff_terms(qpoly)]; alpha=exps[0]
 eqs=[tuple(e[i]-alpha[i] for i in range(n)) for e in exps[1:]]
 rows=[];rhs=[]
 for ell,P in normalized:
  d=s.degree(P,v[p]);d=int(d)
  poly=s.Poly(P,v[p]); R=s.expand(sum(poly.nth(k)*A**(d-k)*(U-B)**k for k in range(d+1)))
  pp=s.Poly(R,U)
  for jj in range(pp.degree()+1):
   Rj=pp.nth(jj)
   if Rj==0:continue
   for ex,c0 in coeff_terms(Rj):
    row=tuple(ell[i]+ex[i]-d*aexp[i]+jj*alpha[i] for i in range(n))
    pair=(row,-jj)
    if pair not in list(zip(rows,rhs)):rows.append(row);rhs.append(-jj)
 cert=solve(rows,rhs,eqs)
 records.append({'factor':str(qpoly),'factor_terms':[[list(e),str(c)] for e,c in k], 'pivot':p+1,'rows':[list(z) for z in rows], 'rhs':rhs,'homogeneity':[list(z) for z in eqs],**cert})
 if idx%10==0: print('factor',idx+1,'bound',cert['upper_bound'],flush=True)
output={'functions':[str(f) for f in func],'nfunctions':len(func),'subset_checks':subsets,'nfactors':len(factors),'zero_case':{'rows':[list(z) for z in allrows],**zero},'cases':records,'factorizations':factorizations}
path=Path(__file__).with_name('e7_divisor_certificates.json');path.write_text(json.dumps(output,indent=2))
print('PASS', len(records),'cases','bounds', {u:sum(r['upper_bound']==u for r in records) for u in ['0','1']},'max denominator',max(r['denominator_max'] for r in records),'bytes',path.stat().st_size,flush=True)