"""Check T494's stored values and the moved recurrence specialization.

Run: cat "$NUMBERDB_KEY_FILE" | agents/sage.sh check_t494_repair.py
Or: NUMBERDB_API_KEY=... sage -python check_t494_repair.py
Reads the live table; never writes it. Outputs reproducible JSON evidence.
"""
import json
import os
import sys
import urllib.request
from decimal import Decimal
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_mpfr import RealField
from sage.version import version as sage_version
import mpmath as mp

key = os.environ.get('NUMBERDB_API_KEY') or sys.stdin.read().strip()
req = urllib.request.Request('https://numberdb.org/api/table?id=T494',
                             headers={'Authorization': 'Bearer ' + key})
with urllib.request.urlopen(req, timeout=60) as response:
    doc = json.load(response)
C = ComplexBallField(600)
R = RealField(600)
mp.mp.dps = 170

def L(nu, x):
    nu, x = C(nu), C(x)
    a = C(QQ(3)/2)
    value = ((x/2)**(nu+1) * (x*x/4).hypergeometric([1], [a, nu+a])
             / (a.gamma() * (nu+a).gamma()))
    assert value.imag().contains_zero()
    assert all(abs(R(part.center())) < R('1e100') and R(part.rad()) < R('1e100')
               for part in (value.real(), value.imag()))
    real = value.real()
    assert real.rad() < R('1e-160') * max(R(1), abs(real.center()))
    return real

def decimal_interval(s):
    d=Decimal(s)
    t=d.as_tuple()
    coeff=QQ(int(''.join(map(str,t.digits)))) * (-1 if t.sign else 1)
    ulp=QQ(10)**int(t.exponent)
    return coeff*ulp,ulp,len(t.digits)

numbers=doc['Numbers']
print(json.dumps({'phase':'shape','orders':list(numbers),'first_group_type':type(next(iter(numbers.values()))).__name__}),flush=True)
rows=[]
for ns,group in numbers.items():
    nu=QQ(ns)
    for xs,entry in group.items():
        s=entry['number'] if isinstance(entry,dict) else entry
        rows.append((nu,QQ(xs),s))
assert len(rows)==600
arguments=sorted(set(x for nu,x,s in rows))
assert len(arguments)==120
expected=sorted(set(QQ(a)/b for b in range(1,5) for a in range(1,20*b+1)))
assert arguments==expected
max_ulp=mp.mpf(0)
for nu,x,s in rows:
    mid,ulp,digits=decimal_interval(s)
    assert digits==100
    ball=L(nu,x)
    assert ball > C(mid-ulp).real() and ball < C(mid+ulp).real()
    mn=mp.mpf(int(nu.numerator()))/int(nu.denominator())
    mx=mp.mpf(int(x.numerator()))/int(x.denominator())
    delta=abs(mp.struvel(mn,mx)-mp.mpf(s))/mp.mpf(str(Decimal(1).scaleb(Decimal(s).as_tuple().exponent)))
    assert delta<1
    max_ulp=max(max_ulp,delta)
print(json.dumps({'phase':'all_values','checked':len(rows),'max_mpmath_discrepancy_ulps':mp.nstr(max_ulp,30)}),flush=True)
max_radius=R(0)
for x in arguments:
    residual=L(QQ(-1),x)-L(QQ(1),x)-2/C.pi().real()
    assert residual.contains_zero()
    assert abs(R(residual.center())) < R('1e100')
    assert residual.rad()<R('1e-160')
    max_radius=max(max_radius,residual.rad())
    # A deliberately false constant must be rejected by the same comparison.
    assert not (L(QQ(-1),x)-L(QQ(1),x)-3/C.pi().real()).contains_zero()
    # The moved formula must also agree with the stored order-one interval.
    mid,ulp,_=decimal_interval(next(s for nu,xx,s in rows if nu==1 and xx==x))
    derived=L(QQ(-1),x)-2/C.pi().real()
    assert derived>C(mid-ulp).real() and derived<C(mid+ulp).real()
print(json.dumps({'sage_version':sage_version,'mpmath_version':mp.__version__,
 'working_bits':600,'mpmath_dps':170,'stored_rows_checked':len(rows),
 'arguments_checked':len(arguments),'negative_order_checks':len(arguments),
 'negative_controls_rejected':len(arguments),
 'negative_order_residual_radius_max':str(max_radius),
 'max_mpmath_discrepancy_ulps':mp.nstr(max_ulp,40),
 'derivation':'In the absolutely convergent defining series for L_-1, the k=0 term is 2/pi. Reindexing k=j+1 makes every remaining term exactly the j-th term of L_1, since the two gamma factors are exchanged.',
 'all_checks_passed':True}),flush=True)
