"""Independent repair checks for numberdb.org/T507, all twelve stored values.

Run beside the attached generate.py with SageMath and numberdb installed:
    sage -python repair_checks.py
For a private draft set NUMBERDB_API_KEY, or pipe the key to stdin with
NUMBERDB_KEY_FROM_STDIN=1. This script only reads tables and computes.

At a=pi*sqrt(r), sum the 16 theta-two terms n=0,...,15. The first omitted
term is 2*exp(-a*(16+1/2)^2); later ratios are at most exp(-34*a), giving
an explicit geometric upper bound. Theta-four uses n=1,...,16 with the
absolute omitted tail bounded by 2*exp(-289*a)/(1-exp(-35*a)).
"""
import json
import math
import os
import sys
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from sage.env import SAGE_VERSION
from generate import ThetaTwoCM

if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
    numberdb.configure(api_key=sys.stdin.read().strip())
B = RealBallField(512)

def exact(s):
    f=Fraction(str(s));return QQ(f.numerator)/QQ(f.denominator)

def interval(s):
    mantissa, _, exponent=s.lower().partition('e')
    places=len(mantissa.partition('.')[2])
    unit=QQ(10)**(int(exponent or 0)-places)
    center=exact(s)
    return center,unit

def finite(x):
    return all(math.isfinite(float(y)) for y in (x.mid(),x.rad()))

def contained(x,s):
    center,unit=interval(s)
    return finite(x) and bool(x>B(center-unit)) and bool(x<B(center+unit))

def value(entry):
    return entry['number'] if isinstance(entry,dict) else entry

def theta2(r):
    a=B.pi()*B(r).sqrt()
    s=sum((2*(-a*(QQ(2*n+1)**2/4)).exp() for n in range(16)),B(0))
    tail=2*(-a*(QQ(33)**2/4)).exp()/(1-(-34*a).exp())
    return s.add_error(tail.upper())

def theta4(r):
    a=B.pi()*B(r).sqrt()
    s=B(1)+sum((2*(-1)**n*(-a*n*n).exp() for n in range(1,17)),B(0))
    tail=2*(-289*a).exp()/(1-(-35*a).exp())
    return s.add_error(tail.upper())

assert interval('0.125')==(QQ(1)/8,QQ(1)/1000)
assert contained(B(QQ(1)/8),'0.125')
assert not contained(B(QQ(1)/8)+QQ(1)/100,'0.125')
assert not contained(B('nan'),'0.125')

doc=numberdb.table('T507')
target=numberdb.table('Jacobi_theta₄_constants_at_singular_nomes')
assert list(doc['Numbers'])==[str(r) for r in range(1,13)]
generator=ThetaTwoCM()
assert generator.digits==100
results=[]
for r,entry in doc['Numbers'].items():
    s=value(entry);x=theta2(int(r));generated=generator.value({'r':int(r)},100)
    y=generated['number'] if isinstance(generated,dict) else generated
    assert contained(x,s) and contained(y,s)
    assert finite(y) and x.overlaps(B(y))
    center,unit=interval(s)
    assert not contained(x+3*B(unit),s)
    assert B(x.rad())<B('1e-150')
    results.append({'r':int(r),'error_ulps_upper':str((abs(x-B(center))/B(unit)).upper()),
                    'radius':str(x.rad()),'generator_accuracy_bits':int(y.accuracy())})
x=theta2(1);y=theta4(1);g=B(QQ(1)/4).gamma()/(2*B.pi())**(QQ(3)/4)
assert finite(y) and finite(g) and x.overlaps(y) and x.overlaps(g)
assert contained(y,value(doc['Numbers']['1']))
assert contained(x,value(target['Numbers']['1']))
assert value(doc['Numbers']['1'])==value(target['Numbers']['1'])
print(json.dumps({'table':'T507','sage':SAGE_VERSION,'working_bits':512,
    'checked_entries':len(results),'unchanged_generator_digits':generator.digits,
    'finite_containment_and_shifted_controls_passed':True,
    'square_point_theta2_theta4_gamma_passed':True,'entries':results},indent=2))
report = generator.verify(sample=None)
print(report)
assert report.ok