"""External high-precision and companion-table checks for numberdb.org/T455."""
import json
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.rational_field import QQ
import mpmath as mp
from checks import endpoints
w=Path('/work')
rows=json.loads((w/'entries.json').read_text())
lookup={(int(e['params']['d1']),int(e['params']['d2']),e['params']['p']):e['number'] for e in rows}
mp.mp.dps=150
p=mp.mpf(19)/20
cdf=lambda x:mp.betainc(4,mp.mpf(1)/2,0,8*x/(8*x+1),regularized=True)
x=mp.findroot(lambda x:cdf(x)-p,(238,239))
lo,hi=endpoints(lookup[(8,1,'19/20')])
def mpq(q): return mp.mpf(int(q.numerator()))/int(q.denominator())
assert mpq(lo)<x<mpq(hi)
print('NIST F(8,1), p=19/20: mpmath 150-digit inversion confirms',mp.nstr(x,115))
assert abs(x-mp.mpf('238.882'))>mp.mpf('0.0005')
assert abs(x-mp.mpf('238.883'))<mp.mpf('0.0005')
# Test the exact Programs snippet by assigning its final expression.
doc=json.loads((w/'document.json').read_text());code=doc['Programs']['program-python']['code']
lines=code.splitlines();lines[-1]='result = '+lines[-1]
namespace={};exec('\n'.join(lines),namespace)
lo,hi=endpoints(lookup[(5,10,'19/20')])
assert mpq(lo)<namespace['result']<mpq(hi)
print('mpmath Programs snippet agrees with the stored-format interval.')
# Companion-table quantiles were generated independently by the t-table run.
# Read its saved database document, keeping the source fixture with the check.
trows=json.loads((w/'T454.json').read_text())['Numbers']
tlookup={(int(e['params']['nu']),e['params']['p']):e['number'] for e in trows}
R=RealBallField(700);I=RealIntervalField(700);count=0
for (d1,d2,p),value in lookup.items():
 if d1!=1:continue
 key=(d2,str((1+QQ(p))/2))
 if key not in tlookup:continue
 x=R(I(*endpoints(value)));t=R(I(*endpoints(tlookup[key])))
 assert x.is_finite() and t.is_finite()
 delta=x-t*t
 assert delta.contains_zero() and delta.abs()<R(10)**(-95)
 count+=1
print('Independent T454 quantiles squared:',count,'overlapping 100-digit intervals.')
assert count==30