"""Checks for the Struve function table -- numberdb.org/T492.

Run with SageMath and numberdb installed:
    sage -python check_t492_repair.py < /path/to/numberdb-key
For a public table stdin may be empty. A saved T492-document.json beside
this script is used if present. This script writes no API data. Its JSON
output records independent decimal comparisons and the negative-order
substitution used to derive the table's Bessel-integral formula.

DLMF 10.22.2 at nu=0 gives pi*x/2*(J0*H_minus1-J_minus1*H0).
Since J_minus1=-J1 and H_minus1=2/pi-H1, substitution gives the
formula printed in T492. The antiderivative tends to zero at the origin.
The Struve recurrence at nu=0 gives H_minus1+H1=1/(sqrt(pi)*Gamma(3/2));
Gamma(3/2)=sqrt(pi)/2 reduces that constant to 2/pi.
"""
import sys
import json
import pathlib
import urllib.request
from decimal import Decimal
import importlib.metadata
import math
import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
import mpmath as mp

p = pathlib.Path(__file__).parent / 'T492-document.json'
if p.exists():
    doc = json.loads(p.read_text())
else:
    key = sys.stdin.read().strip()
    req = urllib.request.Request('https://numberdb.org/api/table?id=T492',
        headers={'Authorization': 'Bearer ' + key} if key else {})
    with urllib.request.urlopen(req, timeout=60) as response:
        doc = json.load(response)

mp.mp.dps = 170
C = ComplexBallField(560)
def mpq(s):
    q = QQ(s)
    return mp.mpf(str(q.numerator())) / mp.mpf(str(q.denominator()))
def struve_ball(nu, x):
    nu, x, a = C(QQ(nu)), C(QQ(x)), C(QQ(3)/2)
    v = (x/2)**(nu+1) * (-(x*x)/4).hypergeometric([1], [a,nu+a]) / (a.gamma()*(nu+a).gamma())
    assert math.isfinite(float(v.real().mid())) and math.isfinite(float(v.imag().mid()))
    assert v.imag().contains_zero()
    r = v.real()
    assert math.isfinite(float(r.rad())) and float(r.rad()) < 1e-140, 'enclosure too wide'
    return r

assert abs(mp.quad(lambda t:t,[0,1])-mp.mpf('0.5')) < mp.mpf('1e-160')
assert abs(mp.quad(lambda t:mp.sin(t),[0,mp.pi/2])-1) < mp.mpf('1e-160')
count=0
worst=mp.mpf(0)
for nu, rows in doc['Numbers'].items():
    for x, written in rows.items():
        assert isinstance(written,str)
        ulp=mp.power(10,int(Decimal(written).as_tuple().exponent))
        err=abs(mp.struveh(mpq(nu),mpq(x))-mp.mpf(written))/ulp
        assert err < 1, (nu,x,str(err))
        worst=max(worst,err); count+=1

constant=1/(C.pi().sqrt()*C(QQ(3)/2).gamma())
assert (constant-2/C.pi()).contains_zero()
args=list(doc['Numbers']['0'])
max_radius=0
for x in args:
    hneg=struve_ball(-1,x); hpos=struve_ball(1,x)
    residual=hneg+hpos-(2/C.pi()).real()
    assert math.isfinite(float(residual.mid())) and residual.contains_zero()
    assert math.isfinite(float(residual.rad())) and float(residual.rad()) < 1e-140
    max_radius=max(max_radius,float(residual.rad()))

max_integral_error=mp.mpf(0)
selected=['1/4','1','13/3','5','10','20']
for s in selected:
    x=mpq(s)
    rhs=mp.pi*x/2*(mp.besselj(0,x)*(2/mp.pi-mp.struveh(1,x))+mp.besselj(1,x)*mp.struveh(0,x))
    intermediate=mp.pi*x/2*(mp.besselj(0,x)*mp.struveh(-1,x)-mp.besselj(-1,x)*mp.struveh(0,x))
    integral=mp.quad(lambda t:mp.besselj(0,t),[0,x/4,x/2,3*x/4,x])
    err=max(abs(rhs-integral),abs(rhs-intermediate))
    assert err < mp.mpf('1e-155'), (s,str(err))
    max_integral_error=max(max_integral_error,err)

print(json.dumps({
    'tid':'T492', 'all_checks_passed':True,
    'versions':{n:importlib.metadata.version(n) for n in ['numberdb','mpmath']},
    'sage_version': __import__('sage.version',fromlist=['version']).version,
    'mpmath_decimal_precision':mp.mp.dps, 'ball_precision_bits':560,
    'stored_decimal_intervals_checked':count,
    'maximum_error_in_last_written_place':mp.nstr(worst,20),
    'negative_order_identity_ball_checks':len(args),
    'maximum_identity_residual_radius':max_radius,
    'bessel_integral_arguments':selected,
    'maximum_bessel_integral_absolute_error':mp.nstr(max_integral_error,12),
    'controls':['integral t on [0,1] = 1/2','integral sin(t) on [0,pi/2] = 1'],
    'source':'https://dlmf.nist.gov/10.22#E2',
    'certification_note':'Mpmath interval comparisons are agreement checks; ball checks enclose the negative-order identity residual.'
},indent=2))