"""Independent checks for numberdb.org/T658.

Run: sage -python checks.py (beside numerics.jsonl)
No generator code is imported. The full-range comparison integrates a direct
quadratic-root Jensen reduction of P_k as a polynomial in y, using x=e^(it).
The generator instead integrates a real trigonometric expression over a
half-period. Algebraic identities and PARI elliptic L-values supply checks
outside this shared Jensen principle. Exterior hypergeometric evaluations
use Rogers--Zudilin (7) only on its stated domain k outside [-4,2].
"""
import sys,json,time
from pathlib import Path
import numberdb.sage as numberdb
import mpmath as mp
from sage.libs.pari import pari
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ

mp.mp.dps=130
HERE=Path(__file__).resolve().parent
rows=[json.loads(x) for x in (HERE/'numerics.jsonl').read_text().splitlines()]
V={r['k']:mp.mpf(r['at190']) for r in rows if 'k' in r}

def argument(key):
    if 'sqrt' not in key: return mp.mpf(key)
    a={'4+2*sqrt(5)':(4,2),'4-2*sqrt(5)':(4,-2),
       '8+4*sqrt(5)':(8,4),'8-4*sqrt(5)':(8,-4)}[key]
    return a[0]+a[1]*mp.sqrt(5)

def compare(label,actual,expected,digits=110):
    assert mp.isfinite(actual) and mp.isfinite(expected)
    error=abs(actual-expected)/max(1,abs(expected))
    assert error<mp.power(10,-digits),(label,mp.nstr(error,8))
    print(json.dumps({'check':label,'relative_error':mp.nstr(error,8),'required_digits':digits}),flush=True)

# A known, nonzero control before the elliptic computations.
pari.set_real_precision(140)
lz=mp.mpf(str(pari('lfun(1,2)')))
compare('PARI lfun control zeta(2)',lz,mp.pi**2/6,120)
L3=(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
mu=3*mp.sqrt(3)*L3/(4*mp.pi)
compare('factorization k=-1',V['-1'],2*mu,120)
compare('spanning tree k=8',V['8'],5*mu,120)
assert V['0']==0

R=PolynomialRing(QQ,names=('x','y'));x,y=R.gens()
P=(1+x)*(1+y)*(x+y)
assert P+x*y==(x+y+1)*(x*y+x+y)
assert P==(1+x)*(1+y)*(x+y)
print(json.dumps({'check':'exact factorizations k=0,-1','passed':True}),flush=True)

# Three independent elliptic L-values, with conductors checked rather than assumed.
curves={14:[1,0,1,4,-6],20:[0,0,0,-108,297],36:[0,0,0,0,1]}
Ls={}
for N,ainvs in curves.items():
    E=pari(ainvs).ellinit()
    assert int(E.ellglobalred()[0])==N
    Ls[N]=mp.mpf(str(pari.lfun(E,2,precision=470)))
    print(json.dumps({'check':'elliptic L value','conductor':N,'ainvs':ainvs,'L2':mp.nstr(Ls[N],125)}),flush=True)
for key,c in [('1',mp.mpf(7)/2),('7',21),('-8',35)]:
    compare('conductor 14 k='+key,V[key],c*Ls[14]/mp.pi**2,120)
compare('conductor 20 k=4',V['4'],10*Ls[20]/mp.pi**2,120)
compare('conductor 20 k=-2',V['-2'],15*Ls[20]/mp.pi**2,120)
compare('conductor 36 k=2',V['2'],9*Ls[36]/(2*mp.pi**2),120)
compare('positive quadratic identity',2*V['4+2*sqrt(5)']-V['8+4*sqrt(5)'],V['4'],120)
assert abs(V['1']-V['-1'])>mp.mpf('0.1')
# This control must fail: a perturbation larger than a stored unit is detected.
assert abs(V['4']+mp.mpf('1e-90')-10*Ls[20]/mp.pi**2)>mp.mpf('1e-100')
print(json.dumps({'check':'deliberate 1e-90 perturbation rejected','passed':True}),flush=True)

# Direct quadratic root integrator; no scalar expression from generate.py.
def roots_integral(key):
    k=argument(key)
    ends=[mp.mpf(0),mp.pi]
    if -1<=k<0:
        ends.extend(2*mp.acos((1+s*mp.sqrt(1+k))/2) for s in [-1,1])
    if 0<k<8:
        ends.append(2*mp.acos((mp.sqrt(1+k)-1)/2))
    ends=sorted(set(ends))
    def integrand(t):
        x=mp.exp(mp.j*t)
        a=1+x; b=(1+x)**2-k*x; c=x*(1+x)
        disc=mp.sqrt(b*b-4*a*c)
        r1=(-b+disc)/(2*a);r2=(-b-disc)/(2*a)
        return mp.log(abs(a))+mp.log(max(1,abs(r1)))+mp.log(max(1,abs(r2)))
    return mp.quad(integrand,ends,method='tanh-sinh',maxdegree=9)/mp.pi
with mp.workdps(125):
    compare('quadratic integrator nonzero control k=-1',roots_integral('-1'),2*mu,115)
    for key in V:
        t=time.monotonic()
        result=roots_integral(key)
        compare('direct quadratic Jensen k='+key,result,V[key],110)

# An independently implemented hypergeometric continuation; real part is used
# only on the source's permitted domain. Boundary k=8 uses its factorization
# control instead of a slowly convergent hypergeometric series at z=1.
with mp.workdps(125):
    def H(z):
        if abs(z) < 1:
            term=mp.mpf(1); total=term; n=0
            while True:
                term *= z*(n+mp.mpf(4)/3)*(n+mp.mpf(5)/3)*(n+1)/(n+2)**3
                total += term
                n += 1
                if abs(term)/(1-abs(z)) < mp.mpf('1e-123'):return total
                if n > 2000000:raise ArithmeticError('Hypergeometric series too slow')
        return mp.hyper([mp.mpf(4)/3,mp.mpf(5)/3,1,1],[2,2,2],z)
    for key in V:
        k=argument(key)
        if not (k < -4 or k > 2) or k==8:continue
        result=mp.re(mp.log((k+4)*(k-2)**4/k**2)-2*k*k/(k+4)**3*H(27*k*k/(k+4)**3)-8*k/(k-2)**3*H(27*k/(k-2)**3))/3
        compare('RZ equation 7 k='+key,result,V[key],110)
    n32=mp.log(32)/3-H(mp.mpf(27)/32)/16
    compare('Hessian normalization g(4)=3*N(32)/4',V['4'],3*n32/4,110)
print(json.dumps({'check':'all independent checks passed','entries':len(V)}),flush=True)