"""Check T522's Euler-gamma notation and unchanged three-step values.

Run with SageMath and numberdb from a directory containing T522-before.json
and skill-hashes.json:
    sage -python repair-gamma-check.py
No network or table writes. Gamma arguments 1/3 and 2/3 are exact rationals;
Arb encloses their conversion and propagates uncertainty throughout.
"""
import json
from decimal import Decimal
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField

here = Path(__file__).resolve().parent
doc = json.loads((here / 'T522-before.json').read_text())
instructions = json.loads((here / 'skill-hashes.json').read_text())

def decimal_parts(text):
    t = Decimal(text).as_tuple()
    coefficient = int(''.join(map(str, t.digits))) * (-1 if t.sign else 1)
    unit = QQ(10) ** t.exponent
    return coefficient * unit, unit

def narrow(value):
    assert value.is_finite() and value > 0
    assert value.rad() < abs(value.mid()) / QQ(10) ** 120
    return value

results = []
for bits in (448, 640):
    R = RealBallField(bits)
    C = ComplexBallField(bits)
    third = QQ(1) / 3
    pi = R.pi()
    # Known controls before any comparison of stored moments.
    assert (narrow(R(1).gamma()) - 1).contains_zero()
    half_gamma_squared = narrow(R(QQ(1) / 2).gamma() ** 2)
    assert half_gamma_squared.overlaps(pi)
    product = narrow(R(third).gamma() * R(2 * third).gamma())
    reflection = narrow(2 * pi / R(3).sqrt())
    assert product.overlaps(reflection)
    A = narrow(3 * R(2) ** third * R(third).gamma() ** 6 / (16 * pi ** 4))
    B = narrow(A + 27 * R(2) ** (2 * third) * R(2 * third).gamma() ** 6 / (4 * pi ** 4))
    # Independent special-function route: elliptic rather than gamma seeds.
    m = (2 - R(3).sqrt()) / 4
    Kc = C(m).elliptic_k()
    assert Kc.real().is_finite() and Kc.imag().is_finite()
    assert Kc.imag().contains_zero()
    K = narrow(Kc.real())
    elliptic_A = narrow(2 * R(3).sqrt() * K ** 2 / pi ** 2)
    elliptic_B = narrow(elliptic_A + R(3).sqrt() / K ** 2)
    assert A.overlaps(elliptic_A) and B.overlaps(elliptic_B)
    moments = {-1: A, 1: B}
    for s in (-1, 1, 3):
        moments[s + 4] = narrow((2 * (5*s*s + 30*s + 46) * moments[s+2]
                                - 9 * (s+2)**2 * moments[s]) / (s+4)**2)
    checked = []
    for order, entry in doc['Numbers']['3'].items():
        text = entry['number'] if isinstance(entry, dict) else entry
        center, ulp = decimal_parts(text)
        value = narrow(moments[int(order)])
        assert abs(value - R(center)) < R(ulp)
        checked.append({'n': 3, 's': int(order), 'enclosed_in_stored_decimal': True,
                        'absolute_error_in_last_places_upper': str((abs(value-R(center))/R(ulp)).upper()),
                        'ball_relative_radius_upper': str((R(value.rad())/abs(value)).upper())})
    results.append({'working_bits': bits, 'controls_passed': True,
                    'gamma_and_elliptic_seeds_overlap': True, 'entries': checked})
print(json.dumps({'table': 'T522', 'instruction_manifest': instructions['instruction_manifest'],
                  'instruction_manifest_sha256': instructions['instruction_manifest_sha256'],
                  'skill_versions': instructions['versions'],
                  'source': 'https://arminstraub.com/downloads/pub/walks2.pdf, equations (4), (5), (7)',
                  'new_values': 0, 'checked_scope': 'Five existing three-step values; four-step values unchanged and not recomputed in this repair.',
                  'results': results}, indent=2))