"""Verify the Freudenthal initial Laurent chart and first-star inclusions.

Requires SymPy. All computations use rational coefficients and all 18
frozen coordinates remain independent. The data are the Laurent expressions
in the auxiliary seed (h,b,c,d,e,f,g,r,q,u) printed in the appendix.
This does not prove that every linear coordinate is a cluster variable.
"""
from pathlib import Path
import json
import sympy as sp
from sympy.polys.fields import field
from check_freudenthal_quadrics import RELATIONS, shift

HERE=Path(__file__).resolve().parent
NAMES=[f'a{i}' for i in range(18)]+list('hbcdefg')+['r','q','u']
F,*xx=field(' '.join(NAMES),sp.QQ)
z=dict(zip(NAMES,xx))
a=[z[f'a{i}'] for i in range(18)]
h,b,c,d,e,f,g,r,q,u=(z[v] for v in list('hbcdefg')+['r','q','u'])
locals_={name:sp.Symbol(name) for name in NAMES}
coords={v:F.from_expr(sp.sympify(value,locals=locals_)) for v,value in json.loads((HERE/'freudenthal_source_coordinates.json').read_text()).items()}
assert len(coords)==56
# Exact defining equations, including all cyclic translates.
for relation in RELATIONS:
    for k in range(18):
        residual=sum((co*coords[shift(x,k)]*coords[shift(y,k)] for co,x,y in relation),F.zero)
        assert not residual, (relation,k)
print('PASS: all 144 cyclic quadratic relations.',flush=True)
# Initial coordinates and the three global higher-degree functions.
for label,var in zip(['x1','x11','x9','y10','x6','x14','x4'],[h,b,c,d,e,f,g]):
    assert coords[label]==var
assert coords['x17']*d-a[17]*c==r
assert e*coords['y13']-a[6]*f==q
assert f*g*coords['y7']-a[14]*g*e-a[3]*f*c-a[4]*f*d==u
for v,value in coords.items():
    assert len(value.denom)==1,(v,'nonmonomial denominator')
    denexp=next(iter(value.denom))
    assert all(denexp[i]==0 for i in range(18)),(v,'frozen denominator')
print('PASS: inverse chart coordinates and polynomial frozen coefficients.',flush=True)
P=a[6]*a[11]*f+a[7]*b*f+a[9]*a[14]*e+a[10]*e*f+a[12]*d*e+a[13]*c*e
Fc0=f*g*(a[6]*a[11]+a[7]*b+a[8]*d+a[10]*e)+a[9]*(a[4]*d*f+a[14]*e*g)
Fr0=g*(a[16]*d+a[17]*b)*h*P+b*q*(d*g*(a[0]*f+a[1]*a[14])+a[17]*f*(a[3]*c+a[4]*d))+h*d*g*q*(a[11]*a[16]+a[15]*b)
ex={
 'h':e*f*q*(a[0]*d*g+a[3]*r+a[17]*(a[3]*c+a[4]*d))+d*f*(e*r+u)*(a[1]*a[6]+a[2]*e)+a[1]*d*q*(u+a[14]*e*g),
 'b':a[11]*(a[6]*f+q)+e*(a[9]*a[14]+a[10]*f+a[12]*d+a[13]*c),
 'c':Fc0+a[9]*u,
 'd':r+a[17]*c,
 'e':q+a[6]*f,
 'f':u+a[14]*e*g,
 'g':u+a[3]*c*f+a[4]*d*f,
 'r':Fr0+g*(a[16]*d+a[17]*b)*b*u+b*q*a[17]*u,
 'q':h*P+b*(e*r+u),
 'u':c*Fr0+b*r*(c*f*(d*(a[1]*a[6]+a[2]*e)+a[3]*q)+h*(Fc0+a[5]*c*d*f+a[13]*c*e*g))
}
norm=dict(zip(list('hbcdefg')+['r','q','u'],[d*e*f,e,f,F.one,F.one,F.one,F.one,b*g,F.one,b*d]))
v=coords
partners={
 'h':v['x2']*q-a[7]*(a[2]*v['x14']+a[1]*v['y13'])-(a[1]*a[6]+a[2]*v['x6'])*v['x12'],
 'b':v['x12'], 'c':v['x4']*v['x8']-a[3]*a[9],
 'd':v['x17'],'e':v['y13'],
 'f':v['x4']*v['y7']-a[4]*v['y10']-a[3]*v['x9'],
 'g':v['x14']*v['y7']-a[14]*v['x6'],
 'r':v['x14']*(v['y2']*v['y7']-v['x3']*v['y10']-v['x1']*v['x8']+a[3]*v['x11'])-v['x6']*(a[13]*v['x1']+a[14]*v['y2']),
 'q':v['x11']*v['y0']-a[11]*v['x1'],
 'u':v['x9']*(v['x1']*v['y16']-a[2]*v['x14']-a[1]*v['y13'])-a[9]*v['x1']*v['x17']
}
for name in ex:
    assert z[name]*norm[name]*partners[name]==ex[name],name
print('PASS: all ten first-mutation identities with their normalizers.',flush=True)
checks=[]
for label,value in coords.items():
    dexp=next(iter(value.denom))
    for name,exchange in ex.items():
        power=dexp[NAMES.index(name)]
        if not power:continue
        assert exchange.denom==F.ring.one
        # The exchange is independent of its own variable. Test each coefficient
        # slice; mere divisibility of the whole numerator would be stronger and wrong.
        slices={}
        i=NAMES.index(name)
        for mon,co in value.numer.items():
            exp=list(mon);j=exp[i]-power;exp[i]=0
            if j<0:
                sl=slices.setdefault(j,F.ring.zero.copy())
                sl[tuple(exp)]=sl.get(tuple(exp),sp.QQ.zero)+co
        for j,sl in slices.items():
            quotient,remainder=sl.div(exchange.numer**(-j))
            assert not remainder,(label,name,j)
            checks.append({'coordinate':label,'direction':name,'pole_order':-j})
print('PASS: all adjacent-chart coefficient divisibilities.',len(checks),flush=True)
report={'status':'PASS','coordinates':56,'frozen_parameters':18,'initial_mutables':10,
 'quadratic_identities':144,'inverse_chart_checks':10,'first_mutation_identities':10,
 'adjacent_negative_slices':len(checks),'max_pole_order':max(c['pole_order'] for c in checks),
 'scope':'Initial Laurent chart, polynomial frozen coefficients, defining quadrics, and first-star membership. Not ordinary-LPA generation or source prime-factor tests.'}
(HERE/'freudenthal_source_report.json').write_text(json.dumps(report,indent=2)+'\n')
(HERE/'freudenthal_source_divisibilities.json').write_text(json.dumps(checks,indent=2)+'\n')
print(json.dumps(report,indent=2))