"""Independent numerical checks for numberdb.org/T630.

Run: sage -python checks.py
For database readback: NUMBERDB_CHECK_STORED=1 sage -python checks.py
An unpublished draft requires NUMBERDB_API_KEY, or NUMBERDB_KEY_FROM_STDIN=1.

The generator uses a polynomial recurrence and an elementary angular-momentum
formula. This checker uses explicit factorial coefficients, validated quadrature
of the defining oblate integrals, the surface equipotential condition, and volume
integrals of mass and axial inertia. Nonzero controls run before comparisons.
The checker borrows only the certified root brackets and candidate values from
the generator. It independently certifies the original equation's endpoint signs.
"""
import json
import os
import platform
from functools import lru_cache
from math import factorial
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 enclose, polynomials, client_from_stdin

BITS = 700
C = ComplexBallField(BITS)
F = RealBallField(BITS)
COARSE_J = ('0.021741','0.029152','0.034638','0.039037','0.042740',
            '0.045957','0.048814','0.051395','0.053755','0.055935')
COARSE_E = ('0.98523','0.99375','0.99657','0.99784','0.99851',
            '0.99891','0.99917','0.99935','0.99947','0.99957')
COARSE_W = ('0.087262','0.066105','0.052711','0.043714','0.037301',
            '0.032513','0.028806','0.025854','0.023449','0.021451')
MANIFEST = 'agents/runs/20261030T195647Z-build-skills/manifest.json'
SKILL_SHA = 'dedecbdca15a728b6e594bfa79f0d149896fca336080ac27b53f513e0b16cd2e'


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


def written_interval(text):
    assert 'e' not in text.lower() and '.' in text
    center = decimal_rational(text)
    ulp = QQ(1)/10**len(text.split('.')[1])
    return center-ulp, center+ulp


@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):
    v=x.parent()(0)
    for a in explicit_coefficients(n)[::-1]:
        v=v*x*x+a
    return v


def integral_q(n,x,threshold=150):
    # t=x/u maps [x,infinity] to [0,1]. No special-function recurrence.
    coeff=[a*x**(2*j) for j,a in enumerate(explicit_coefficients(n))]
    def integrand(u,analytic):
        p=C(0)
        for a in coeff:
            p=p*u*u+a
        return x*u**(2*n)/((u*u+x*x)*p*p)
    v=explicit_p(n,x)*C.integral(integrand,0,1)
    assert v.real().is_finite() and v.imag().is_finite()
    assert v.imag().contains_zero()
    assert F(v.real().rad()) < F(10)**(-threshold)
    return v.real()


def original_equation(k,x):
    return (2*explicit_p(2*k+2,x)*integral_q(2*k+2,C(x))
            -x/explicit_p(2,x)*(1-2*x*integral_q(2,C(x))))


def independent_physics(x,scale=1,rho=1,gravity=1):
    # Focal length is scale; equatorial squared radius a2 and polar radius c.
    scale,rho,gravity=F(scale),F(rho),F(gravity)
    a2=scale**2*(1+x*x)
    c=scale*x
    q0=integral_q(0,C(x),threshold=115)
    q2=integral_q(2,C(x),threshold=115)
    pref=-F(4)/3*F.pi()*gravity*rho*scale**2*x*(1+x*x)
    phi_p=pref*(q0+q2)
    phi_e=pref*(q0-q2/2)
    omega_squared=2*(phi_e-phi_p)/a2
    # M = pi rho a^2 c int_-1^1 (1-z^2) dz.
    mass_integral=C.integral(lambda z,analytic: 1-z*z,-1,1).real()
    inertia_integral=C.integral(lambda z,analytic: (1-z*z)**2,-1,1).real()
    mass=F.pi()*rho*a2*c*mass_integral
    inertia=F.pi()/2*rho*a2*a2*c*inertia_integral
    assert mass > 0 and inertia > 0 and omega_squared > 0
    j2=inertia**2*omega_squared*rho**(QQ(1)/3)/(4*F.pi()*gravity*mass**(QQ(10)/3))
    assert j2.is_finite() and F(j2.rad()) < F(10)**(-115)
    return j2, omega_squared/(4*F.pi()*gravity*rho), inertia/(mass*a2)


def contained(v,lo,hi):
    assert v.is_finite() and v > F(lo) and v < F(hi)


def main(document=None):
    # Known nonzero controls test the quadrature before any bifurcation values.
    for n,target in [(0,F.pi()/4),(2,F.pi()/2-F(3)/2)]:
        v=integral_q(n,C(1))
        assert v.overlaps(target) and not v.contains_zero()
    for power,target in [(1,QQ(4)/3),(2,QQ(16)/15)]:
        v=C.integral(lambda z,analytic:(1-z*z)**power,-1,1).real()
        assert v.is_finite() and F(v.rad())<F(10)**(-150)
        assert v.overlaps(F(target)) and not v.contains_zero()
    print('CONTROLS PASS: q0(1), q2(1), mass and inertia integrals.',flush=True)
    stored=os.environ.get('NUMBERDB_CHECK_STORED')=='1' or document is not None
    if stored:
        doc=document or numberdb.table('T630',client=client_from_stdin())
        records=doc['Numbers']
        if isinstance(records,dict):
            assert list(map(str,records))==list(map(str,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]==list(map(str,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':'T630','stored_readback':stored,'skill_manifest':MANIFEST,
              'local_skill_sha256':SKILL_SHA,'live_skill_sha256':SKILL_SHA,
              'sage_version':SAGE_VERSION,'python_version':platform.python_version(),
              'working_bits':500,'independent_bits':BITS,'controls_passed':True,'entries':[]}
    previous=F(0)
    for k in range(1,11):
        j,lo,hi=enclose(k)
        j700,lo700,hi700=enclose(k,working_bits=700)
        assert (lo,hi)==(lo700,hi700) and j.overlaps(j700)
        p,_,_,_,num=polynomials(k)
        assert tuple(p.list()[::2])==explicit_coefficients(2*k+2)
        left=original_equation(k,F(lo));right=original_equation(k,F(hi))
        assert left>0 and right<0
        x=F((lo+hi)/2).add_error(F((hi-lo)/2))
        independent,omega2,inertia_ratio=independent_physics(x)
        rescaled,_,_=independent_physics(x,scale=2,rho=3,gravity=5)
        assert independent.overlaps(rescaled) and independent.overlaps(j)
        assert inertia_ratio.overlaps(F(QQ(2)/5))
        # Every displayed conversion, with propagated root uncertainty.
        e=1/(1+x*x).sqrt();s=(1-e*e).sqrt();theta=e.arcsin()
        formula_omega=(s*(3-2*e*e)*theta/(e**3)-3*s*s/(e*e))/2
        formula_j=QQ(4)/25*(3/(4*F.pi()))**(QQ(4)/3)*formula_omega/(1-e*e)**(QQ(2)/3)
        beta=3/(2*e*e)-1-3*s/(2*e*theta)
        beta_j=QQ(4)/25*(3/(4*F.pi()))**(QQ(4)/3)*beta*theta/(e*s**(QQ(1)/3))
        for v in [formula_omega,formula_j,beta_j]:
            assert v.is_finite() and F(v.rad())<F(10)**(-115)
        assert omega2.overlaps(formula_omega)
        assert j.overlaps(formula_j) and j.overlaps(beta_j)
        a,b=written_interval(texts[k])
        assert len(texts[k].split('.')[1])==101
        for v in [j,j700,independent,rescaled,formula_j,beta_j]:contained(v,a,b)
        assert previous<j
        previous=j
        for v,text in [(j,COARSE_J[k-1]),(e,COARSE_E[k-1]),(omega2,COARSE_W[k-1])]:
            a0,b0=written_interval(text);mid=(a0+b0)/2;half=(b0-a0)/4
            contained(v,mid-half,mid+half)
        row={'k':k,'j2':texts[k],'radius':str(j.rad()),'independent_radius':str(independent.rad()),
             'xi_lower':str(lo),'xi_upper':str(hi),'original_integral_endpoint_residuals':[str(left),str(right)],
             'uniqueness_polynomial_degree':int(num.degree()),'same_bracket_at_700_bits':True,
             'written_interval_contains_all_checks':True,'scale_density_G_invariance':True,
             'all_displayed_conversions':True,'coarse_source_checks':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()