"""Recheck T493's unchanged numbers during its presentation repair.

Run with SageMath and a JSON table document obtained from GET /api/table?id=T493:
    sage -pip install numberdb
    sage -python repair_checks.py < T493-document.json

This is a read-only check. It regenerates all 360 values at 600 bits, checks
containment in the stored decimal intervals, and independently tests the
integral at three arguments using elementary ball quadrature at 480 bits.
It does not replace the table's original generator or its full certification.
"""
import json
import sys
from decimal import Decimal
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION


def decimal_interval(text):
    value = Decimal(text)
    midpoint = QQ(Fraction(value))
    ulp = QQ(10) ** int(value.as_tuple().exponent)
    assert len(value.as_tuple().digits) == 100, (text, value.as_tuple())
    return midpoint, ulp


def regenerated_values(doc, bits=600):
    C = ComplexBallField(bits)
    i = C(0, 1)
    rotation = (1 + i) / C(2).sqrt()
    for order, rows in doc['Numbers'].items():
        n = int(order)
        for argument, text in rows.items():
            value = ((-i)**n * (C(QQ(argument))*rotation).bessel_K(n)).real()
            yield order, argument, text, value


def main():
    doc = json.load(sys.stdin)
    assert list(doc['Parameters']) == ['nu', 'x']
    assert list(doc['Numbers']) == ['0', '1', '2']
    expected = sorted({QQ(a)/b for b in range(1, 5) for a in range(1, 20*b+1)})
    assert len(expected) == 120
    for rows in doc['Numbers'].values():
        assert [QQ(x) for x in rows] == expected
    max_radius = 0
    max_ulp_error = 0
    count = 0
    for n, x, text, ball in regenerated_values(doc):
        midpoint, ulp = decimal_interval(text)
        R = ball.parent()
        assert ball.is_finite(), (n, x, 'non-finite')
        assert R(ball.rad()) < R(ulp)/1000000, (n, x, 'wide ball')
        error = abs(ball - R(midpoint)) / R(ulp)
        assert error.upper() < 1, (n, x, 'outside written interval')
        max_ulp_error = max(max_ulp_error, float(error.upper()))
        max_radius = max(max_radius, float(ball.rad()))
        count += 1
    assert count == 360
    print(json.dumps({'regenerated_rows': count, 'precision_bits': 600,
                      'max_absolute_radius': max_radius,
                      'max_error_in_written_ulps': max_ulp_error,
                      'all_rows_have_100_significant_digits': True}), flush=True)

    C = ComplexBallField(480)
    R = RealBallField(480)
    tolerance = R(2)**(-440)
    controls = []
    for f, start, end, expected in [
        (lambda t, analytic: t, C(0), C(1), QQ(1)/2),
        (lambda t, analytic: t.sin(), C(0), C.pi()/2, QQ(1)),
    ]:
        answer = C.integral(f, start, end, abs_tol=tolerance, rel_tol=tolerance)
        assert answer.real().is_finite() and answer.imag().is_finite() and (answer.real()-R(expected)).contains_zero() and answer.imag().contains_zero()
        assert R(answer.real().rad()) < R(2)**(-400)
        controls.append(str(answer))

    b = R(1)/(8*R(2).sqrt())
    u0 = R(9).exp()
    tail = (-b*u0).exp()*(u0/b + 1/b**2)
    assert tail.upper() < R(QQ(819)/100 * QQ(10)**(-307)).lower()
    integral_checks = []
    for n, x in [(0, QQ(1)/4), (1, QQ(1)), (2, QQ(20))]:
        a = C(x)/C(2).sqrt()
        phase = n*C.pi()/2
        def integrand(t, analytic):
            z = a*t.cosh()
            return (-z).exp()*(z+phase).cos()*(n*t).cosh()
        value = C.integral(integrand, 0, 9, abs_tol=tolerance, rel_tol=tolerance)
        assert value.real().is_finite() and value.imag().is_finite() and value.imag().contains_zero()
        midpoint, ulp = decimal_interval(doc['Numbers'][str(n)][str(x)])
        error = abs(value.real()-R(midpoint)) + tail
        assert R(value.real().rad()) < R(ulp)/1000000
        assert error.upper() < R(ulp).lower(), (n,x,'integral mismatch')
        integral_checks.append({'nu':n, 'x':str(x), 'radius':str(value.real().rad()),
                                'contained_with_tail':True})

    i = C(0,1)
    half = QQ(1)/2
    value = ((-i*C.pi()*C(half)/2).exp()*((1+i)/C(2).sqrt()).bessel_K(half)).real()
    elementary = (R.pi()/2).sqrt()*(-1/R(2).sqrt()).exp()*(1/R(2).sqrt()+3*R.pi()/8).cos()
    assert value.is_finite() and elementary.is_finite()
    assert R(value.rad()) < R(2)**(-400) and R(elementary.rad()) < R(2)**(-400)
    assert value.overlaps(elementary)
    print(json.dumps({'integral_checks':integral_checks, 'control_integrals':controls,
                      'tail_upper':str(tail.upper()), 'half_order_control':True,
                      'sage_version':SAGE_VERSION,
                      'numberdb_version':getattr(numberdb,'__version__','unknown')}),flush=True)


if __name__ == '__main__':
    main()