"""Numerical checks for T541, separate full-trace formulas and controls.

Run: sage -python checks.py [--stored]
Download generate.py, accelerate.py, and check.py alongside this file.
The high-precision Voros check shares only the generic summation helper with
the Watkins computation. An optional --library check instead uses mpmath's
own hypergeometric implementation at 25 decimal digits on every N>=3.
Neither high-precision agreement nor that library check proves a tail bound.
"""
import os,sys,json,time
import numberdb.sage as numberdb
import mpmath as mp
from sage.rings.real_mpfi import RealIntervalField
from generate import compute,boundary,FullSpectralZeta
from accelerate import hyper_tail
from check import agrees_with,measure,exactness,names_its_rings

# This formula is transcribed independently from Voros 1983 (C.22)-(C.28).
def voros(N,n,working=190,library=False):
    with mp.workdps(working):
        if N==1:
            # Independently take log derivatives of the Airy determinants.
            return mp.nstr(-mp.diff(lambda E:mp.log(mp.airyai(-E)*(-mp.airyai(-E,1))),0,2),working)
        if N==2:
            # Spectral odd-integer Hurwitz-zeta representation.
            return mp.nstr(mp.zeta(n,mp.mpf(1)/2)/mp.mpf(2)**n,working)
        v=mp.mpf(1)/(N+2);h=mp.mpf(1)/2
        parity=mp.sin(mp.pi*v)/(2*mp.sqrt(mp.pi))*(2*v)**(2-4*v)
        parity*=mp.gamma(v)*mp.gamma(2*v)*mp.gamma(3*v)/mp.gamma(2*v+h)
        if n==1:
            return mp.nstr(mp.tan(2*mp.pi*v)/mp.tan(mp.pi*v)*parity,working)
        A=[4*v,2*v,3*v,5*v,h]
        B=[1+2*v,1+v,1-v,h+4*v]
        f=mp.hyper(A,B,1) if library else hyper_tail(A,B,240,160)
        pref=mp.sin(mp.pi*v)/(4*v*v*mp.sqrt(mp.pi))*(2*v)**(4-8*v)
        pref*=mp.gamma(3*v)*mp.gamma(4*v)*mp.gamma(5*v)/mp.gamma(h+4*v)
        return mp.nstr(parity**2+pref*f,working)


def close(a,b,tolerance):
    return abs(a-b)<tolerance*max(mp.mpf(1),abs(b))


def run_checks(stored_mode=False,library=False):
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    mp.mp.dps=220
    gen=FullSpectralZeta()
    # Known, nonzero control precedes any interpretation of the method.
    with mp.workdps(190):
        control=hyper_tail([mp.mpf(1)]*2,[mp.mpf(3)],240,160)
        assert abs(control-2)<mp.mpf('1e-150')
    assert not close(mp.mpf(2)+mp.mpf('1e-90'),mp.mpf(2),mp.mpf('1e-110'))
    print('CONTROL 2F1(1,1;3;1)=2 and deliberate 1e-90 perturbation: PASS',flush=True)
    from numberdb._generate import _stored_entries
    records=_stored_entries(numberdb.table('T541')) if stored_mode else None
    field=RealIntervalField(800)
    vals={};independent={};comparison={};reference={};minimum=mp.inf
    for p in gen.enumerate():
        N,n=p['N'],p['n'];key=f'N={N},n={n}'
        actual=mp.mpf(compute(N,n,200))
        other=mp.mpf(voros(N,n,25 if library else 190,library))
        agreement=-mp.log10(abs((actual-other)/actual)) if actual!=other else mp.inf
        assert agreement>(20 if library else 115),(key,'Voros disagreement',str(agreement))
        assert actual>0
        independent[key]=field(str(other))
        vals[key]=gen.value(p,100)
        compare_digits=18 if library else 100
        comparison[key]=mp.nstr(actual,compare_digits)
        reference[key]=mp.nstr(other,compare_digits)
        if N>=3:
            # Watkins 2.10a,b, derived from the quantum Wronskian, checks both
            # parity sectors independently of the generic summation formula.
            z1m=mp.mpf(boundary(N,1,1,200));z1p=mp.mpf(boundary(N,1,-1,200))
            sig=mp.mpf(2)/(N+2);lam=mp.mpf(1)/2
            coeff=lambda a:mp.sin(mp.pi*sig*(lam+a))/mp.sin(mp.pi*sig*lam)
            if n==1:
                residual=coeff(1)*z1m+coeff(-1)*z1p
            else:
                zm=mp.mpf(boundary(N,2,1,200));zp=mp.mpf(boundary(N,2,-1,200))
                residual=coeff(2)*zm+coeff(-2)*zp+(coeff(1)**2-coeff(2))*(z1p-z1m)**2
            assert abs(residual)<mp.mpf('1e-120'),(key,'Wronskian',str(residual))
        if stored_mode and not library:
            text=records[f'{N},{n}']
            from decimal import Decimal
            center=Decimal(text); ulp=Decimal(10)**center.as_tuple().exponent
            mid=field(str(center));rad=field(str(ulp))
            interval=field(mid.lower()-rad.upper(),mid.upper()+rad.upper())
            assert interval.lower()<independent[key].lower() and independent[key].upper()<interval.upper(),key
        minimum=min(minimum,agreement)
        print(json.dumps({'key':key,'independent_agreement_digits':str(agreement),'stored':stored_mode,'library':library, 'wronskian_residual':str(residual) if N>=3 else None}),flush=True)
    assert len(vals)==22
    assert not exactness(vals)
    assert not names_its_rings(__file__.replace('checks.py','generate.py'))
    assert not agrees_with(comparison, lambda key: reference[key])
    print('MINIMUM_INDEPENDENT_AGREEMENT',mp.nstr(minimum,12),'MEASURE',measure(vals),flush=True)
    return vals

if __name__=='__main__':
    run_checks('--stored' in sys.argv,'--library' in sys.argv)