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

Run: sage -python checks.py
Set NUMBERDB_CHECK_STORED=1 to read the draft/public table instead of local
values (an unpublished draft requires NUMBERDB_API_KEY). This script shares
only the generated values with generate.py: p_n is computed from factorial
coefficients, q_n from its defining integral using certified Arb quadrature.
The original AKM equation is evaluated at both endpoints of every written
100-digit eccentricity. The two endpoint signs must differ.
"""
import json
import os
from math import factorial
import platform
from functools import lru_cache
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION
from generate import MaclaurinEccentricities, enclose, polynomials, client_from_stdin

PRECISION = 500
C = ComplexBallField(PRECISION)
F = RealBallField(PRECISION)
COARSE_E = ('0.98523','0.99375','0.99657','0.99784','0.99851',
            '0.99891','0.99917','0.99935','0.99947','0.99957')
COARSE_AXIS = ('0.17126','0.11160','0.082750','0.065744','0.054534',
               '0.046589','0.040664','0.036075','0.032417','0.029432')

def decimal_rational(text):
    if 'e' in text.lower():
        mant, power = text.lower().split('e')
        return decimal_rational(mant)*QQ(10)**int(power)
    a,b = text.split('.')
    return QQ(int(a+b))/QQ(10)**len(b)

@lru_cache(None)
def explicit_coefficients(n):
    assert n % 2 == 0
    return tuple(QQ(factorial(n+2*j))/(2**n*factorial(n//2-j)*
                 factorial(n//2+j)*factorial(2*j)) for j in range(n//2+1))

def explicit_p(n,x):
    value = x.parent()(0)
    for c in explicit_coefficients(n)[::-1]:
        value = value*x*x+c
    return value


def integral_q(n,x):
    # Substitute t=x/u in the integral over [x,infinity].
    coefficients = [a*x**(2*j) for j,a in enumerate(explicit_coefficients(n))]
    def integrand(u,analytic):
        p_scaled = C(0)
        for a in coefficients:
            p_scaled = p_scaled*u*u+a
        return x*u**(2*n)/((u*u+x*x)*p_scaled*p_scaled)
    value = explicit_p(n,x)*C.integral(integrand,0,1)
    assert value.real().is_finite() and value.imag().is_finite()
    assert F(value.real().rad()) < F(10)**(-130)
    assert value.imag().contains_zero()
    return value.real()


def original_equation(k,e):
    x = ((1-e*e)/(e*e)).sqrt()
    xc = C(x)
    return (2*explicit_p(2*k+2,x)*integral_q(2*k+2,xc)
            -x/explicit_p(2,x)*(1-2*x*integral_q(2,xc)))


def sign(v):
    assert v.is_finite()
    if v > 0: return 1
    if v < 0: return -1
    raise ArithmeticError('Unresolved sign in independent integral check.')


def main(document=None):
    # Controls must succeed before results are believed.
    controls = [(0,F.pi()/4),(2,F.pi()/2-F(3)/2)]
    for n,known in controls:
        got = integral_q(n,C(1))
        assert got.overlaps(known)
        assert not got.contains_zero()
    print('CONTROL: q0(1)=pi/4 and q2(1)=pi/2-3/2 certified.', flush=True)
    stored = os.environ.get('NUMBERDB_CHECK_STORED') == '1'
    table = document if document is not None else (
        numberdb.table('T627',client=client_from_stdin()) if stored else None)
    if document is not None:
        assert stored, 'An exported document is a stored-value check.'
    if stored:
        records=table['Numbers']
        if isinstance(records, dict):
            assert [str(k) for k in records]==[str(k) for k in range(1,11)]
            texts={int(k):(v['number'] if isinstance(v,dict) else v)
                   for k,v in records.items()}
        else:
            assert [str(r['params']['k']) for r in records]==[str(k) for k in range(1,11)]
            texts={int(r['params']['k']):r['number'] for r in records}
    else:
        texts={k:to_text(enclose(k)[0],100) for k in range(1,11)}
    evidence={'table':'T627','stored_readback':stored,
              'skill_manifest':'agents/runs/20261030T141911Z-build-skills/manifest.json',
              'skill_sha256':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
              'sage_version':SAGE_VERSION,'python_version':platform.python_version(),
              'working_bits':500,'second_working_bits':700,'entries':[]}
    previous = F(0)
    for k in range(1,11):
        e,lo,hi=enclose(k)
        other,lo2,hi2=enclose(k,working_bits=700)
        assert (lo,hi)==(lo2,hi2) and e.overlaps(other)
        assert previous < e
        previous=e
        p,_,_,_,derivative=polynomials(k)
        assert list(explicit_coefficients(2*k+2)) == list(p.list()[::2])
        center=decimal_rational(texts[k]); ulp=QQ(1)/10**100
        # Each eccentricity lies in (0.1,1), so 100 digits means 100 decimals.
        assert len(texts[k].split('.')[1])==100
        assert e > F(center-ulp) and e < F(center+ulp)
        minus=original_equation(k,F(center-ulp))
        plus=original_equation(k,F(center+ulp))
        assert sign(minus)==-1 and sign(plus)==1
        coarse=decimal_rational(COARSE_E[k-1]); coarse_axis=decimal_rational(COARSE_AXIS[k-1])
        assert abs(e-F(coarse)) < F(QQ(1)/200000)
        axis=(1-e*e).sqrt()
        halfunit=QQ(1)/(2*10**len(COARSE_AXIS[k-1].split('.')[1]))
        assert abs(axis-F(coarse_axis)) < F(halfunit)
        row={'k':k,'eccentricity':texts[k], 'xi_lower':str(lo),'xi_upper':str(hi),
             'eccentricity_radius':str(e.rad()),'independent_endpoint_residuals':[str(minus),str(plus)],
             'global_uniqueness_certificate_degree':int(derivative.degree()),
             'positive_coefficient_pattern_verified':True,
             'coarse_eccentricity_and_axis_ratio_match':True,'same_bracket_at_700_bits':True}
        evidence['entries'].append(row)
        print('PASS',k,texts[k],flush=True)
    print('EVIDENCE_JSON='+json.dumps(evidence,sort_keys=True),flush=True)

if __name__=='__main__':
    main()