"""Independent checks for Even period polynomials -- numberdb.org/T695.
Run: sage -python checks.py. Needs generate.py and optional T694-roster.json.
The Sage exact generator is compared with PARI's independent modular-form
construction and numerical mfperiodpol, and with explicit split Mellin sums.
Numerical agreement is a check, not the source of the exact integer entries.
"""
import json,math,hashlib
from pathlib import Path
from functools import lru_cache
import numberdb.sage as numberdb
from sage.libs.pari import pari
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.real_mpfr import RealField
from sage.rings.complex_mpfr import ComplexField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.qqbar import AA
from sage.arith.misc import euler_phi
from sage.env import SAGE_VERSION
from generate import EvenPeriodPolynomials,space_data,selected_spaces,exact_polynomial,QX,X
from check import exactness,measure,names_its_rings
BASE=Path(__file__).resolve().parent
P=PolynomialRing(QQ,'x')
def pol(c):return P(str(pari.lift(c)).replace('y','x'))
@lru_cache(None)
def dimensions(N,k):
    primes=ZZ(N).prime_divisors()
    mu=QQ(N)
    for p in primes:mu*=1+QQ(1)/p
    e2=0 if N%4==0 else math.prod(1+int(ZZ(-4).kronecker(p)) for p in primes)
    e3=0 if N%2==0 or N%9==0 else math.prod(1+int(ZZ(-3).kronecker(p)) for p in primes)
    cusps=sum(euler_phi(math.gcd(d,N//d)) for d in ZZ(N).divisors())
    genus=1+mu/12-QQ(e2)/4-QQ(e3)/3-QQ(cusps)/2
    full=(k-1)*(genus-1)+(k//2-1)*cusps+e2*(k//4)+e3*(k//3)
    new=full-sum(len(ZZ(N//d).divisors())*dimensions(int(d),k)[1] for d in ZZ(N).divisors() if d<N)
    return int(full),int(new)

def pari_roster(N,k,B):
    mf=pari.mfinit([N,k],0);rows=[]
    for f in pari.mfeigenbasis(mf):
        fp=pol(pari.mfparams(f)[3]);cs=[pol(c) for c in pari.mfcoefs(f,B)]
        for z in fp.roots(AA,multiplicities=False):
            rows.append((tuple(c(z) for c in cs),f,int(fp.degree())))
    rows.sort(key=lambda row:row[0][2:])
    return mf,rows

def eta_delta(nmax):
    # q product_(m>=1) (1-q^m)^24, independent integer convolution.
    a=[1]+[0]*(nmax-1)
    for m in range(1,nmax):
        b=[0]*nmax
        for j in range(0,min(24,(nmax-1)//m)+1):
            c=(-1)**j*math.comb(24,j)
            for n in range(nmax-j*m):b[n+j*m]+=c*a[n]
        a=b
    return [0]+a

def mellin_coefficients(coeff,N,k,epsilon,working,terms):
    R=RealField(numberdb.bits(working,losing=32));T=1/R(N).sqrt();twopi=2*R.pi()
    integrals=[R(0) for s in range(k)]
    for n in range(1,terms+1):
        a=R(coeff[n]);t=twopi*n;z=t*T
        decay=(-z).exp();partial=R(1);term=R(1)
        for s in range(1,k):
            if s>1:term*=z/(s-1);partial+=term
            integrals[s]+=a*decay*math.factorial(s-1)*partial/t**s
    moments={s:integrals[s]+epsilon*R(N)**(k//2-s)*integrals[k-s] for s in range(1,k)}
    w=k-2
    c=[R(0)]*(w+1)
    for j in range(0,w+1,2):c[w-j]=math.comb(w,j)*(-1)**(j//2)*moments[j+1]
    assert c[w]!=0
    return [v/c[w] for v in c]

def exact_identities(p,N,k,fricke):
    w=k-2
    assert p.degree()==w and p.leading_coefficient()>0
    assert math.gcd(*(int(c) for c in p.list()))==1
    assert p(-X)==p(X)
    reciprocal=sum(QQ(p[w-j])*QQ(N)**(j-w//2)*X**j for j in range(w+1))
    assert QX(p)==-fricke*reciprocal
    if N==1:
        manin=p+sum(p[j]*(X-1)**j*X**(w-j) for j in range(w+1))+sum(p[j]*(-1)**j*(X-1)**(w-j) for j in range(w+1))
        assert manin==0 and p(1)==p(-1)==0


def run(values=None):
    # Nonzero controls run first.
    assert eta_delta(5)==[0,1,-24,252,-1472,4830]
    control=list(pari.mfcoefs(pari.mfDelta(),100))
    assert [int(c) for c in control]==eta_delta(100)
    delta=36*X**10-691*X**8+2073*X**6-2073*X**4+691*X**2-36
    assert exact_polynomial(1,12,1)==delta
    R=RealField(numberdb.bits(150));T=R(1);a=2*R.pi()
    assert abs((-a*T).exp()/a-R('0.000297212741'))<R('1e-12')
    # Explicit antiderivative control for nonzero exponential integral.
    assert abs(((-a*T).exp()/a)*a*(-a*T).exp()**(-1)-1)<R('1e-140')
    report={'sage':SAGE_VERSION,'pari':str(pari('version()')),'controls':{'Delta_eta_coefficients':100,'Delta_published_even_polynomial':True,'nonzero_exponential_integral':True},'spaces':[],'entries':[],'companion_indices_checked':0}
    vals={} if values is None else values
    companion_path=BASE/'T694-roster.json'
    companion=json.loads(companion_path.read_text()) if companion_path.exists() else []
    comp={(r['N'],r['k'],r['i']):r for r in companion}
    for N,k in selected_spaces():
        M,C,S,rows=space_data(N,k);B=max(2,int(M.sturm_bound()))
        full,new=dimensions(N,k)
        assert (full,new)==(C.dimension(),S.dimension())
        mf,other=pari_roster(N,k,B)
        assert int(pari.mfdim(mf))==new and int(pari.mfdim([N,k],1))==full
        assert len(other)==len(rows)
        for r,pr in zip(rows,other):
            assert r['sort']==pr[0] and r['field_degree']==pr[2]
            key=(N,k,r['i'])
            if comp:
                old=comp[key];fp=P(old['field_poly']);lo,hi=map(QQ,old['root_interval'])
                zs=[z for z in fp.roots(AA,multiplicities=False) if lo<z<hi];assert len(zs)==1
                cc=tuple(P([QQ(c) for c in cs])(zs[0]) for cs in old['coefficients'])
                assert cc==r['sort'][:len(cc)] and old['field_degree']==r['field_degree']
                report['companion_indices_checked']+=1
        report['spaces'].append({'N':N,'k':k,'cusp_dimension':full,'new_dimension':new,'rational_forms':sum(r['field_degree']==1 for r in rows)})
        for r,pr in zip(rows,other):
            if r['field_degree']!=1:continue
            i=r['i'];key=(N,k,i);f=pr[1];w=k-2
            p=exact_polynomial(*key) if values is None else QX(vals[key])
            vals[key]=p
            Bcheck=max(100,B)
            qs=r['factor'].q_eigenform(Bcheck+1,names='a')
            pc=[QQ(c) for c in pari.mfcoefs(f,Bcheck)]
            assert list(qs.padded_list(Bcheck+1))==pc
            atkin=pari.mfatkininit(mf,N,precision=128);assert atkin[2]==1
            ac=list(pari.mfcoefs(pari.mfatkin(atkin,f),B))
            fc=list(pari.mfcoefs(f,B))
            fricke=next(s for s in [-1,1] if all(a==s*b for a,b in zip(ac,fc)))
            exact_identities(p,N,k,fricke)
            epsilon=(-1)**(k//2)*fricke
            analytic=[]
            for working in [160,220]:
                bits=numberdb.bits(working,losing=32);pari.set_real_precision(working+12);CC=ComplexField(bits);RR=CC.base_ring()
                raw=pari.mfperiodpol(mf,f,1,precision=bits)
                coef=[CC(pari.polcoef(raw,j)) for j in range(w+1)]
                assert abs(coef[w])>RR('1e-20')
                norm=[c/coef[w] for c in coef]
                err=max(abs(norm[j]-CC(QQ(p[j])/p[w])) for j in range(w+1))
                assert not err.is_NaN() and not err.is_infinity() and err<RR(10)**(-working+30),(key,working,err)
                analytic.append({'working_digits':working,'maximum_absolute_coefficient_ratio_error':str(err)})
            coeff=[QQ(c) for c in pari.mfcoefs(f,600)]
            mellin=[]
            for working,terms in [(150,450),(190,600)]:
                mc=mellin_coefficients(coeff,N,k,epsilon,working,terms)
                RR=mc[0].parent()
                err=max(abs(mc[j]-RR(QQ(p[j])/p[w])) for j in range(w+1))
                assert not err.is_NaN() and not err.is_infinity() and err<RR(10)**(-working+25),(key,working,err)
                mellin.append({'working_digits':working,'terms':terms,'maximum_absolute_coefficient_ratio_error':str(err)})
            report['entries'].append({'N':N,'k':k,'i':i,'fricke':fricke,'epsilon':epsilon,'pari':analytic,'mellin':mellin,'exact_identities':True,'q_coefficients_crosschecked':Bcheck})
            print('CHECKED',*key,'PARI',analytic[-1]['maximum_absolute_coefficient_ratio_error'],'Mellin',mellin[-1]['maximum_absolute_coefficient_ratio_error'],flush=True)
    assert len(vals)==18 and not exactness(vals)
    assert not names_its_rings(BASE/'generate.py')
    report['measure']=measure(vals)
    report['manifest']='agents/runs/20261109T174608Z-build-skills/manifest.json'
    report['manifest_sha256']=hashlib.sha256((BASE/'manifest.json').read_bytes()).hexdigest()
    report['skill_sha256']='886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b'
    print('EVIDENCE '+json.dumps(report))
    return report
if __name__=='__main__':run()