"""Independent checks of T691. Run: sage -python check_values.py.

Requires generate.py and numerical-evidence.json in the same directory.
Uses Sage modular-symbol eigenforms, the square-divisor coefficient formula,
and Dokchitser's evaluator independently of the generator's exact Eisenstein
series and reciprocal Euler polynomials. Also checks Rankin unfolding at the
right endpoint against PARI's Haberland Petersson pairing for every embedding.
Reports actual differences; neither evaluator certifies an analytic remainder.
"""
import os
import json
from pathlib import Path
from math import isqrt
import numberdb.sage as numberdb
import sage.symbolic.expression
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
from sage.rings.complex_mpfr import ComplexField
from sage.modular.modform.constructor import Newforms
from sage.lfunctions.dokchitser import Dokchitser
from sage.libs.pari import pari
from generate import exact_forms, tensor_coefficients, objects

BITS=700
R=RealField(BITS)
C=ComplexField(BITS)


def relative(a,b):
    return abs(a-b)/abs(b)


def run(evidence=None):
    if evidence is None:
        evidence=json.loads(Path(__file__).with_name('numerical-evidence.json').read_text())
    records={tuple(int(r['params'][n]) for n in ('k','i','s')):r for r in evidence['records']}
    report={'dokchitser_bits':BITS,'controls':{},'objects':[],'comparisons':[]}
    # Nonzero numerical controls precede all measurements.
    z=Dokchitser(conductor=1,gammaV=[0],weight=1,eps=1,poles=[1],residues=[-1],init='1',prec=BITS)
    zv=z(2); zerr=relative(zv,C(R.pi()**2/R(6)))
    assert zv.real()>1 and zerr<R(10)**-180
    report['controls']['dokchitser_zeta2_relative_error']=str(zerr)
    assert list(pari.mfcoefs(pari.mfDelta(),4))==[0,1,-24,252,-1472]
    report['controls']['pari_delta_coefficients']=[0,1,-24,252,-1472]
    # Confirm roster dimensions independently of the proposed list.
    dims={k:int(pari.mfdim([1,k],0)) for k in range(2,28,2)}
    assert {k:d for k,d in dims.items() if d}=={12:1,16:1,18:1,20:1,22:1,24:2,26:1}
    report['dimensions']=dims
    for k,i in objects():
        D=Dokchitser(conductor=1,gammaV=[0,1,-11,-10],weight=k+11,eps=1,prec=BITS)
        n=int(D.num_coeffs())+30
        F=Newforms(1,k,names='a')[0]
        G=Newforms(1,12,names='d')[0]
        K=F.base_ring()
        embeddings=sorted(K.embeddings(R),key=lambda phi:phi(F[2]))
        phi=embeddings[i-1]
        f=F.q_expansion(n+1); g=G.q_expansion(n+1)
        forms=exact_forms(n)
        # Exact equality includes the Sturm range and all numerical inputs.
        assert all(f[m]==K(forms[k][m][0])+K(forms[k][m][1])*F[2] for m in range(n+1))
        assert all(g[m]==forms[12][m][0] for m in range(n+1))
        for p in [2,3,5,7]:
            for m in range(1,k//12+1):
                assert f[p*m]+(ZZ(p)**(k-1)*f[m//p] if m%p==0 else 0)==f[p]*f[m]
        coeff=[K(0)]*(n+1)
        for d in range(1,isqrt(n)+1):
            d2=d*d; power=ZZ(d)**(k+10)
            for m in range(1,n//d2+1): coeff[d2*m]+=power*f[m]*g[m]
        primary=tensor_coefficients(k,n)
        assert all(coeff[m]==K(primary[m][0])+K(primary[m][1])*F[2] for m in range(1,n+1))
        # Inputs are formed at more precision than the requested Dokchitser result.
        RR=RealField(BITS+256)
        hi_phi=sorted(K.embeddings(RR),key=lambda ph:ph(F[2]))[i-1]
        D.init_coeffs([hi_phi(a) for a in coeff[1:]])
        maxrel=R(0); maximag=R(0)
        vals={}
        for s in range(12,k):
            v=D(s); target=R(records[k,i,s]['high']); err=relative(v,C(target))
            assert (not v.real().is_NaN() and not v.real().is_infinity()) and (not v.imag().is_NaN() and not v.imag().is_infinity())
            assert err<R(10)**-130,(k,i,s,err)
            maxrel=max(maxrel,err); maximag=max(maximag,abs(v.imag()))
            vals[s]=v
            report['comparisons'].append({'k':k,'i':i,'s':s,'value':str(v),'relative_error':str(err)})
        feq=D.check_functional_equation()
        # Check the explicit reflected identity on independently computed values.
        worst_fe=R(0)
        for s in range(12,k):
            t=k+11-s
            ratio=(2*R.pi())**(2*(t-s))*ZZ(s-1).factorial()*ZZ(s-12).factorial()/(ZZ(t-1).factorial()*ZZ(t-12).factorial())
            err=relative(vals[t],vals[s]*ratio)
            assert err<R(10)**-130
            worst_fe=max(worst_fe,err)
        # The endpoint oracle uses a period-pairing algorithm, not lfun.
        pari.set_real_precision_bits(BITS)
        mf=pari.mfinit([1,k],0); pf=pari.mfeigenbasis(mf)[0]
        pg=pari.mfmul(pari.mfDelta(),pari.mfEk(k-12))
        fs=pari.mfsymbol(mf,pf,precision=BITS)
        gs=pari.mfsymbol(mf,pg,precision=BITS)
        pairing=pari.mfpetersson(fs,gs)
        if k==24:
            emb=list(pari.mfembed(pf,pari.mfcoefs(pf,2)[2],precision=BITS))
            order=sorted(range(2),key=lambda h:R(emb[h]))
            inner=C(pairing[order[i-1]])
        else:
            inner=C(pairing)
        endpoint=inner*C(pari.zeta(k-12,precision=BITS))*(4*R.pi())**(k-1)/ZZ(k-2).factorial()
        enderr=relative(endpoint,C(R(records[k,i,k-1]['high'])))
        assert abs(inner)>R('1e-8') and enderr<R(10)**-130,(k,i,enderr)
        item={'k':k,'i':i,'a2':str(phi(F[2])),'coefficients':n,
              'exact_coefficients_match':True,'hecke_sturm_checks':True,
              'max_relative_difference':str(maxrel),'max_imaginary_part':str(maximag),
              'dokchitser_functional_equation_error':str(feq),
              'max_reflected_relative_error':str(worst_fe),
              'petersson_pairing':str(inner),'endpoint_from_pairing':str(endpoint),
              'endpoint_relative_difference':str(enderr)}
        report['objects'].append(item)
        print('CHECKED',k,i,'Dokchitser',maxrel,'Petersson',enderr,flush=True)
    print('CHECKS_JSON '+json.dumps(report),flush=True)
    return report


if __name__=='__main__': run()