"""Independent checks for Scorer Gi values -- numberdb.org/T514.

Run: sage -python checks.py
Read-back checks: sage -python checks.py --stored (key required for the draft).

The production generator sums the Gi Taylor series. This checker instead
uses rigorous quadrature of the entire Airy functions in DLMF 9.12.4:
  Gi(x) = Bi(x)*(1/3 - integral_0^x Ai(t)dt) + Ai(x)*integral_0^x Bi(t)dt.
The constant integral_0^infinity Ai(t)dt = 1/3 is DLMF 9.10.11.
Both paths use Arb, but the functions and integration algorithm differ.
The analytic flag may be ignored because Ai and Bi are entire.

Controls precede comparisons: the same quadrature integrates t^2 to 1/3,
and Ai'(t) integrates to Ai(1)-Ai(0). Every compared ball must be finite
and narrow. The mpmath comparisons at 150 and 200 decimal digits are
heuristic independent checks, not the basis for the table's proven label.
The second derivative in the ODE and Hi in Gi+Hi=Bi come from mpmath;
neither is derived from the production Gi values.

Stdout ends with EVIDENCE_JSON followed by machine-readable results.
"""
import json
import os
import sys
import time
import importlib.metadata

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
from numberdb._write import to_text
import mpmath as mp

from generate import ARGUMENTS, scorer_gi, ScorerGiValues
from check import exactness, measure, names_its_rings


def require_narrow(value, bits=380):
    assert value.is_finite(), value
    assert value.accuracy() >= bits, value
    return value


def airy_integral(x):
    C = ComplexBallField(665)
    x = C(QQ(x))
    tol = QQ(1)/2**510
    ia = C.integral(lambda t, analytic: t.airy_ai(), 0, x,
                    abs_tol=tol, rel_tol=tol)
    ib = C.integral(lambda t, analytic: t.airy_bi(), 0, x,
                    abs_tol=tol, rel_tol=tol)
    result = x.airy_bi()*(C(QQ(1)/3)-ia) + x.airy_ai()*ib
    assert result.real().is_finite() and result.imag().is_finite()
    assert result.imag().contains_zero()
    return require_narrow(result.real())


def mp_argument(text):
    q = QQ(text)
    return mp.mpf(int(q.numerator()))/int(q.denominator())


def mp_midpoint(ball):
    q = QQ(ball.center())
    return mp.mpf(int(q.numerator()))/int(q.denominator())


def written_interval(text):
    # The generator writes ordinary decimals. Interpret their final place
    # exactly, without conversion through a binary float or Sage QQ('0.1').
    from decimal import Decimal
    d = Decimal(text)
    numerator, denominator = d.as_integer_ratio()
    centre = QQ(numerator)/denominator
    ulp = QQ(10)**d.as_tuple().exponent
    return centre-ulp, centre+ulp


def contained(ball, text):
    lo, hi = written_interval(text)
    assert ball > lo and ball < hi, (ball, text)


def main():
    started = time.monotonic()
    C = ComplexBallField(665)
    R = RealBallField(665)
    tol = QQ(1)/2**510
    control = C.integral(lambda t, analytic:t*t, 0, 1,
                         abs_tol=tol, rel_tol=tol).real()
    require_narrow(control)
    assert (control-R(QQ(1)/3)).contains_zero()
    airy_control = C.integral(lambda t, analytic:t.airy_ai_prime(), 0, 1,
                              abs_tol=tol, rel_tol=tol).real()
    airy_expected = (C(1).airy_ai()-C(0).airy_ai()).real()
    require_narrow(airy_control)
    require_narrow(airy_expected)
    assert airy_control.overlaps(airy_expected)
    print('Controls passed: polynomial integral and Airy derivative integral.',flush=True)
    evidence = {'table':'T514','controls':{
        'polynomial_accuracy_bits':int(control.accuracy()),
        'airy_derivative_accuracy_bits':int(airy_control.accuracy())},
        'versions':{'python':sys.version,'numberdb':numberdb.__version__,
                    'mpmath':mp.__version__},'entries':{}}
    try:
        from sage.version import version
        evidence['versions']['sage'] = version
    except ImportError:
        pass
    stored = None
    if '--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
            os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
        stored = numberdb.table('T514')['Numbers']
        print('Read database Numbers:',len(stored),flush=True)
    values = {}
    texts = {}
    for x in ARGUMENTS:
        value, diagnostic = scorer_gi(x, diagnostics=True)
        values[x] = value
        text = to_text(value,100)
        texts[x] = text
        integral = airy_integral(x)
        require_narrow(value)
        assert integral.overlaps(value), x
        contained(value,text)
        contained(integral,text)
        low_high = {}
        for precision in (150,200):
            with mp.workdps(precision):
                xx = mp_argument(x)
                independently = mp.scorergi(xx)
                relative = abs((mp_midpoint(value)-independently)/independently)
                assert relative < mp.mpf('1e-130'), (x,relative)
                low_high[str(precision)] = mp.nstr(independently,precision)
                second = mp.diff(mp.scorergi,xx,2)
                ode = second-xx*mp_midpoint(value)+1/mp.pi
                identity = mp_midpoint(value)+mp.scorerhi(xx)-mp.airybi(xx)
                assert abs(ode) < mp.mpf('1e-130'), (x,ode)
                assert abs(identity) < mp.mpf('1e-130'), (x,identity)
                if precision == 200:
                    diagnostic.update({'mpmath_relative_error':mp.nstr(relative,8),
                                       'ode_residual':mp.nstr(ode,8),
                                       'Gi_plus_Hi_minus_Bi':mp.nstr(identity,8)})
        with mp.workdps(210):
            assert abs(mp.mpf(low_high['150'])-mp.mpf(low_high['200'])) < mp.mpf('1e-140')
        if stored is not None:
            # Current API uses a nested parameter mapping, matching table display order.
            row = stored[x]
            stored_text = row['number'] if isinstance(row,dict) else row
            contained(integral,stored_text)
            contained(value,stored_text)
            with mp.workdps(200):
                centre = sum(written_interval(stored_text))/2
                mcentre = mp.mpf(int(centre.numerator()))/int(centre.denominator())
                lo,hi = written_interval(stored_text)
                radius = (hi-lo)/2
                mradius = mp.mpf(int(radius.numerator()))/int(radius.denominator())
                xx = mp_argument(x)
                assert abs(mp.scorergi(xx)-mcentre) < mradius
                assert abs(mp.diff(mp.scorergi,xx,2)-xx*mcentre+1/mp.pi) <= abs(xx)*mradius+mp.mpf('1e-180')
                assert abs(mcentre+mp.scorerhi(xx)-mp.airybi(xx)) <= mradius+mp.mpf('1e-180')
            diagnostic['stored_text'] = stored_text
            diagnostic['stored_ODE_and_sum_identity'] = True
        diagnostic.update({'written':text,'airy_integral_accuracy_bits':int(integral.accuracy()),
                           'airy_integral_ball':str(integral),'mpmath':low_high})
        evidence['entries'][x] = diagnostic
        print('PASS x='+x+', Taylor bits='+str(value.accuracy())+
              ', Airy-integral bits='+str(integral.accuracy()),flush=True)
    with mp.workdps(200):
        assert abs(mp.scorergi(0)-1/(mp.power(3,mp.mpf(7)/6)*mp.gamma(mp.mpf(2)/3))) < mp.mpf('1e-190')
        assert abs(mp.diff(mp.scorergi,0)-1/(mp.power(3,mp.mpf(5)/6)*mp.gamma(mp.mpf(1)/3))) < mp.mpf('1e-190')
    origin = R(C(0).airy_bi().real()/3)
    require_narrow(origin)
    assert origin.overlaps(values['0'])
    assert not exactness(values)
    assert not names_its_rings(__file__.replace('checks.py','generate.py'))
    evidence['measure_balls'] = measure(values)
    evidence['measure_written'] = measure(texts)
    evidence['seconds'] = time.monotonic()-started
    evidence['stored_checked'] = stored is not None
    print('EVIDENCE_JSON',flush=True)
    print(json.dumps(evidence,indent=2),flush=True)


if __name__ == '__main__':
    main()