"""T619 Programs repair checks. Run: sage -python repair_checks.py.

Download program-repaired.py, generate.py, tidal_solver.py, document-before.json,
and policy-evidence.json into the same directory. This script writes no table.
The analytic cases are checked through a separate hypergeometric logarithmic
 derivative, including both stored degrees and additional degrees.
"""
from pathlib import Path
import json,hashlib,platform,contextlib,io
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.version import version as sage_version
import mpmath as mp
from generate import PolytropicLove
root=Path(__file__).resolve().parent
p=root/'program-repaired.py'
namespace={};printed=io.StringIO()
with contextlib.redirect_stdout(printed):exec(compile(p.read_text(),str(p),'exec'),namespace)
love=namespace['love_analytic']
doc=json.loads((root/'document-before.json').read_text())
assert list(PolytropicLove().enumerate())==[{'n':n,'ell':ell} for n,group in doc['Numbers'].items() for ell in group]
try:PolytropicLove().value({'n':'0','ell':'4'},digits=100)
except ValueError as e:old_refusal=str(e)
else:raise AssertionError('Original refusal absent')

def decimal_rational(s):
    from fractions import Fraction
    f=Fraction(s);return QQ(f.numerator)/f.denominator

def independent(ell,dps):
    with mp.workdps(dps):
        # H(x) proportional to x^ell * 0F1(ell+3/2; -x^2/4).
        # Its logarithmic derivative uses the derivative of 0F1.
        b=mp.mpf(ell)+mp.mpf(3)/2
        f=mp.hyp0f1(b,-mp.pi**2/4)
        fp=mp.hyp0f1(b+1,-mp.pi**2/4)
        y=ell-mp.pi**2*fp/(2*b*f)
        k=(ell-y)/(2*(ell+1+y))
        assert mp.isfinite(k) and k>0
        return mp.nstr(k,dps-10)

results=[]
for ell in list(range(2,31))+[50,100]:
    exact=love(0,ell)
    assert exact==QQ(3)/(4*(ell-1))
    # Also check the modified-field exterior relation using y=ell-3.
    y=QQ(ell-3)
    assert exact==(ell-y)/(2*(ell+1+y))
    got=love(1,ell)
    assert got.is_finite() and got>0 and got.accuracy()>=numberdb.bits(100,losing=8)
    a=decimal_rational(independent(ell,180));b=decimal_rational(independent(ell,220))
    assert abs((a-b)/b)<QQ(1)/10**165
    assert b in got,(ell,got,b)
    rec={'ell':ell,'n0_exact':str(exact),'n1_ball':str(got),'accuracy_bits':int(got.accuracy()),'independent_220_digits':str(b)}
    if ell in (2,3):
        assert exact==QQ(doc['Numbers']['0'][str(ell)]['number'])
        stored=doc['Numbers']['1'][str(ell)]['number']
        from decimal import Decimal
        ulp=QQ(10)**Decimal(stored).as_tuple().exponent
        midpoint=decimal_rational(stored)
        # Every point in the computed ball lies in the written +/-1 ulp interval.
        assert got>midpoint-ulp and got<midpoint+ulp
        assert not (got>midpoint+4*ulp and got<midpoint+6*ulp)
        rec['stored_interval_contains_ball']=True
        rec['midpoint_error_ulps']=str(RealBallField(80)(abs(b-midpoint)/ulp))
    results.append(rec)
assert love(0,4)==QQ(1)/4
for n,ell in [('1/2',4),(1,'5/2'),(1,1)]:
    try:love(n,ell)
    except ValueError:pass
    else:raise AssertionError('Invalid domain input accepted')
evidence={'table':'T619','run':'20261030T061717Z','skill_policy':json.loads((root/'policy-evidence.json').read_text()),
          'source':'https://arxiv.org/html/1709.02469v1#S5.SS3','source_equations':[28,33],
          'derivation':'The spherical Bessel recurrence (2ell+1) j_ell(x)=x[j_(ell-1)(x)+j_(ell+1)(x)] simplifies equation (33) to j_(ell+1)(pi)/(2 j_(ell-1)(pi)); the common sqrt(pi/(2x)) factor cancels in the cylindrical Bessel ratio.',
          'original_refusal':old_refusal,'program_sha256':hashlib.sha256(p.read_bytes()).hexdigest(),
          'example_output':printed.getvalue(),'environment':{'sage':sage_version,'python':platform.python_version(),'mpmath':mp.__version__},
          'checks':results,'all_passed':True,'minimum_accuracy_bits':min(r['accuracy_bits'] for r in results),
          'stored_analytic_entries_checked':4,'all_20_identities_match_generator':True,
          'limits':'No nonlinear value or production generator was changed. This repair checks all four analytic stored entries and 31 degrees of each analytic family, not the 16 unchanged nonlinear values. The independent mpmath check uses two precisions, not a rigorous independent enclosure. Arb propagates pi uncertainty in the replacement program.'}
print(json.dumps(evidence,indent=2))