"""Check T618's stored intervals and cross-table identities after filling.
Run: sage -python readback_checks.py (NUMBERDB_API_KEY needed for drafts).
Reads live tables; writes nothing. Fresh generator verification precedes
independent checks on the actual stored decimal intervals, not their display
midpoints alone. The attached radial evidence is a separate computation.
"""
import json
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfi import RealIntervalField
from generate import LaneEmdenMassRadius,INDICES,coefficient,background

I=RealIntervalField(768)
def stored_interval(v):
    if isinstance(v,dict):v=v['number']
    text=str(v);f=Fraction(text);q=QQ(f.numerator)/QQ(f.denominator)
    ulp=QQ(10)**Decimal(text).as_tuple().exponent if '.' in text or 'e' in text.lower() else QQ(0)
    a=I(q-ulp,q+ulp)
    assert a.lower()>0 and a.absolute_diameter()<I('3e-98')
    return a

def overlaps(a,b):
    return not (a.upper()<b.lower() or b.upper()<a.lower())

def main():
    report=LaneEmdenMassRadius().verify(sample=None)
    print('VERIFY',report,flush=True);assert report.ok
    docs={tid:numberdb.table(tid) for tid in ('T618','T613','T614','T615','T617')}
    values={tid:{str(n):stored_interval(v) for n,v in d['Numbers'].items()} for tid,d in docs.items()}
    assert tuple(values['T618'])==INDICES
    independent=json.loads(Path(__file__).with_name('independent-evidence.json').read_text())['values']
    output={}
    for n in INDICES:
        nn=I(QQ(n));N=values['T618'][n];x=values['T613'][n];q=values['T614'][n]
        C=values['T615'][n];W=values['T617'][n]
        high=coefficient(n,210);a,b,_=background(n,210)
        assert N.lower()<=high<=N.upper(),('stored precision',n)
        assert x.lower()<=a<=x.upper() and q.lower()<=b<=q.upper()
        surface=(4*I.pi())**(1/nn)/(nn+1)*q**((1-nn)/nn)*x**((nn-3)/nn)
        pressure=W*(4*I.pi()/(3*C))**(1+1/nn)
        assert overlaps(N,surface),('surface identity',n)
        assert overlaps(N,pressure),('pressure identity',n)
        for scale,density in ((I(1),I(1)),(I(2),I(3))):
            M=4*I.pi()*scale**3*density*q;R=scale*x
            K=4*I.pi()*scale**2*density**(1-1/nn)/(nn+1)
            assert overlaps(K,N*M**((nn-1)/nn)*R**((3-nn)/nn)),('scale identity',n)
        if n=='1':
            assert overlaps(N,2/I.pi())
            assert not overlaps(N+I('1e-80'),2/I.pi())
            K=4*I.pi()/2;R=x
            assert overlaps(R,(I.pi()*K/2).sqrt())
        if n=='3':
            assert overlaps(N,(4*I.pi())**(I(1)/3)/(4*q**(I(2)/3)))
            M=4*I.pi()*q;K=I.pi()
            assert overlaps(M,(K/N)**(I(3)/2))
        independent_N=I(independent[n]['coefficient'])
        err=abs(independent_N/N-1)
        assert err.upper()<I('1e-40')
        output[n]={'stored_interval_contains_fresh_210_digit_value':True,
                   'surface_and_pressure_identities_overlap':True,
                   'physical_scale_identity_at_two_scales':True,
                   'independent_relative_error_bound':str(err.upper()),
                   'stored_interval_width':str(N.absolute_diameter())}
    print('READBACK '+json.dumps(output),flush=True)
    return output
if __name__=='__main__':main()