"""Independent checks for prolate differential eigenvalues -- numberdb.org/T518.

    sage -python checks.py
    sage -python checks.py --stored

The second command reads T518 and checks every stored decimal interval by
opposite signs of an independent endpoint-series residual. A private draft
requires NUMBERDB_API_KEY (or NUMBERDB_KEY_FROM_STDIN=1 and a piped key).
Output is JSON evidence, suitable for saving as evidence.json. No writes.

The endpoint solution is normalized here by f(1)=1, independently of the
unit-norm convention; eigenvalues do not depend on this normalization.
Writing t=1-x and f=sum a_k t^k in the differential equation gives
  2(k+1)^2 a_(k+1) = (k(k+1)+q-chi)a_k - 2q a_(k-1) + q a_(k-2).
We impose f(0)=0 for odd modes or f'(0)=0 for even modes. This does not
construct a Legendre matrix or share its arithmetic. Fixed-precision series
and sign tests are diagnostic, not rigorous remainder bounds.
"""
import json
import os
import sys
import time
import platform
from pathlib import Path
import importlib.metadata
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
from scipy.special import pro_cv
import scipy
import sage.version
import generate as g


def endpoint(chi, q, parity, terms):
    R = chi.parent()
    a, previous, older = R(1), R(0), R(0)
    da, dprevious, dolder = R(0), R(0), R(0)
    total = R(1) if parity else R(0)
    dtotal = R(0)
    for k in range(terms-1):
        factor = k*(k+1)+q-chi
        divisor = 2*(k+1)**2
        new = (factor*a - 2*q*previous + q*older)/divisor
        dnew = (factor*da - a - 2*q*dprevious + q*dolder)/divisor
        weight = 1 if parity else k+1
        total += weight*new
        dtotal += weight*dnew
        older, previous, a = previous, a, new
        dolder, dprevious, da = dprevious, da, dnew
    return total, dtotal


def endpoint_root(S, n, dps, terms):
    R = RealField(numberdb.bits(dps, losing=16))
    q = (R.pi()*R(S)/2)**2
    # Deliberately only a double-precision seed from an outside implementation.
    chi = R(str(pro_cv(0, n, float(R.pi()*S/2))))
    for iteration in range(12):
        residual, derivative = endpoint(chi, q, n%2, terms)
        assert derivative and (not derivative.is_infinity() and not derivative.is_NaN()) and (not residual.is_infinity() and not residual.is_NaN())
        change = residual/derivative
        chi -= change
        if abs(change) < R(10)**(-(dps-30)):
            return chi, iteration+1
    raise ArithmeticError('Endpoint Newton iteration did not converge')


def agreement_digits(a, b):
    R = RealField(900)
    difference = abs(R(a)-R(b))
    if not difference:
        return 250
    return int(-(difference/max(abs(R(a)),R(1))).log10())


def controls():
    # Exact, nonzero controls before interpreting zero residuals.
    assert endpoint(QQ(0), QQ(0), 1, 20)[0] == 1  # P0(0)
    assert endpoint(QQ(2), QQ(0), 0, 20)[0] == -1 # d P1(1-t)/dt
    assert endpoint(QQ(6), QQ(0), 1, 20)[0] == -QQ(1)/2 # P2(0)
    for n in range(g.MODES):
        residual, derivative = endpoint(QQ(n*(n+1)), QQ(0), n%2, 30)
        assert residual == 0 and derivative != 0
        computed = g.bisect_mode(g.field(160)(0), n, 20, 160)
        assert abs(computed-n*(n+1)) < g.field(160)(10)**(-140)
    # Algebraic check of the DLMF shift, in every degree used in the largest matrix.
    for k in range(360):
        diag = QQ((k+1)**2)/((2*k+1)*(2*k+3)) + QQ(k*k)/((2*k-1)*(2*k+1))
        dlmf = -QQ(2*(k*(k+1)-1))/((2*k-1)*(2*k+3))
        assert diag-1 == dlmf
        upper = -QQ((k+1)*(k+2))/((2*k+3)*(2*k+5))
        lower = -QQ((k+1)*(k+2))/((2*k+1)*(2*k+3))
        square = QQ((k+1)**2*(k+2)**2)/((2*k+1)*(2*k+3)**2*(2*k+5))
        assert upper*lower == square
    return {'endpoint_exact_nonzero_controls':3,'zero_bandwidth_modes':12,
            'dlmf_shift_rational_coefficient_checks':360}


def parse_decimal(text):
    """Exact rational midpoint and last decimal unit, without binary floats."""
    mantissa, marker, exponent = text.lower().partition('e')
    exponent = int(exponent) if marker else 0
    integer, dot, fraction = mantissa.partition('.')
    assert dot
    power = exponent-len(fraction)
    unit = QQ(10)**power
    return QQ(int(integer+fraction, 10))*unit, unit


def stored_values():
    tree = numberdb.table('T518')
    values = {}
    def walk(node, keys):
        if len(keys)==3:
            value = node['number'] if isinstance(node,dict) else node
            values[tuple(map(int, keys))] = value
            return
        if isinstance(node,dict) and 'numbers' in node:
            node=node['numbers']
        for key, value in node.items():
            walk(value, keys+[key])
    nums = tree['Numbers']
    if isinstance(nums,list):
        for item in nums:
            p=item['params']
            values[(int(p['m']),int(p['S']),int(p['n']))]=item['number']
    else:
        walk(nums,[])
    assert len(values)==48
    return values


def run(use_stored=False):
    start=time.monotonic()
    report={'table':'T518','controls':controls(),'settings':g.SETTINGS,
            'endpoint_settings':[(220,600),(260,800)],'rows':[],
            'versions':{'python':platform.python_version(),'sage':sage.version.version,
                        'scipy':scipy.__version__}}
    try:
        report['versions']['numberdb']=importlib.metadata.version('numberdb')
    except importlib.metadata.PackageNotFoundError:
        report['versions']['numberdb']='repository client'
    stored=stored_values() if use_stored else None
    for S in g.BANDWIDTHS:
        spectra=[g.spectrum(str(S),*setting) for setting in g.SETTINGS]
        R=g.field(260)
        q=(R.pi()*S/2)**2
        previous=None
        for n in range(g.MODES):
            reference=R(spectra[-1][n])
            settings_agreement=min(agreement_digits(v[n],reference) for v in spectra[:-1])
            assert settings_agreement>=110
            root1, it1=endpoint_root(S,n,220,600)
            root2, it2=endpoint_root(S,n,260,800)
            independent_agreement=min(agreement_digits(reference,root1),agreement_digits(reference,root2))
            assert independent_agreement>=110
            assert n*(n+1)<reference<n*(n+1)+q
            if previous is not None:
                assert reference>previous
            previous=reference
            scipy_error=abs(reference-R(str(pro_cv(0,n,float(R.pi()*S/2)))))/reference
            assert scipy_error<R('1e-12')
            text=stored[(0,S,n)] if stored is not None else g.decimal_string(reference,100)
            midpoint,ulp=parse_decimal(text)
            assert abs(R(midpoint)-root2)<R(ulp)
            # Direct test of the interval asserted by the written decimal.
            # At both series sizes, it brackets a root by opposite residual signs.
            for dps,terms in ((220,600),(260,800)):
                F=RealField(numberdb.bits(dps,losing=16))
                qf=(F.pi()*S/2)**2
                fl=endpoint(F(midpoint-ulp),qf,n%2,terms)[0]
                fh=endpoint(F(midpoint+ulp),qf,n%2,terms)[0]
                assert (not fl.is_infinity() and not fl.is_NaN()) and (not fh.is_infinity() and not fh.is_NaN()) and fl*fh<0
            report['rows'].append({'m':0,'S':S,'n':n,'number':text,
               'precision_truncation_agree_digits':settings_agreement,
               'endpoint_agree_digits':independent_agreement,
               'endpoint_self_agree_digits':agreement_digits(root1,root2),
               'endpoint_iterations':[it1,it2],
               'scipy_relative_error':str(scipy_error.n(30)),
               'decimal_interval_sign_changes':2,
               'rounding_error_ulps':str((abs(R(midpoint)-root2)/R(ulp)).n(30))})
        print('checked S=%d'%S,file=sys.stderr,flush=True)
    report['read_from_database']=use_stored
    report['minimum_precision_truncation_agreement']=min(r['precision_truncation_agree_digits'] for r in report['rows'])
    report['minimum_endpoint_agreement']=min(r['endpoint_agree_digits'] for r in report['rows'])
    report['seconds']=time.monotonic()-start
    return report

if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    print(json.dumps(run('--stored' in sys.argv),indent=2))