"""Independent exact checks of Debye derivative polynomials -- numberdb.org/T568.

Run: sage -python checks.py. To recheck an API document read back from the
website: sage -python checks.py --stored table.json. No API key is needed.

The independent construction differentiates DLMF 10.32.9 under the integral:
K'_nu(nu*z) = -1/2 integral_R cosh(t) exp[-nu*(z*cosh(t)-t)] dt.
At t0=asinh(1/z), set p=1/sqrt(1+z^2), q=sqrt(p), epsilon=nu^(-1/2),
and t=t0+q*epsilon*s. After removing the leading prefactor, the Gaussian
expectation has amplitude cosh(q*epsilon*s)+q^2*sinh(q*epsilon*s).
The nonquadratic exponent has epsilon^j coefficient
-q^(j + 2*(j mod 2))*s^(j+2)/(j+2)!.
Formal exponentiation and exact Gaussian moments yield (-1)^k V_k(p)
at epsilon^(2k). This uses no U recurrence or derivative conversion.
"""
import sys
import json
import pathlib
from math import factorial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from check import exactness, measure, agrees_with, names_its_rings

P = PolynomialRing(QQ, 'p')
p = P.gen()
S = PolynomialRing(QQ, names=('q', 's'))
q, s = S.gens()

def moment(n):
    if n % 2:
        return QQ(0)
    result = QQ(1)
    for j in range(1, n, 2):
        result *= j
    return result

def gaussian(poly):
    out = P(0)
    for (q_power, s_power), coefficient in poly.dict().items():
        if s_power % 2:
            continue
        assert q_power % 2 == 0
        out += coefficient*moment(s_power)*p**(q_power // 2)
    return out

def independent_polynomials(max_k):
    nmax = 2*max_k
    phase = [S(0)] + [-q**(j+2*(j % 2))*s**(j+2)/QQ(factorial(j+2))
                       for j in range(1, nmax+1)]
    exponential = [S(1)]
    for n in range(1, nmax+1):
        exponential.append(sum((j*phase[j]*exponential[n-j]
                                for j in range(1, n+1)), S(0))/QQ(n))
    amplitude = [(q*s)**n*(q**2 if n % 2 else 1)/QQ(factorial(n))
                 for n in range(nmax+1)]
    result = {}
    for n in range(nmax+1):
        term = sum((amplitude[j]*exponential[n-j] for j in range(n+1)), S(0))
        integrated = gaussian(term)
        if n % 2:
            assert integrated == 0
        else:
            result[n // 2] = (-1)**(n // 2)*integrated
    return result

def validate(values):
    # Controls precede the comparisons: nonzero known Gaussian moments,
    # including the value 1, plus one odd moment.
    assert [moment(j) for j in (0, 1, 2, 4, 6)] == [1, 0, 1, 3, 15]
    assert gaussian(S(1)) == 1 and gaussian(s**2) == 1
    assert gaussian(q**2*s**4) == 3*p
    expected = independent_polynomials(max(values))
    assert expected[0] == 1
    assert expected[1] == (7*p**3-9*p)/QQ(24)
    assert not exactness(values)
    assert not agrees_with(values, lambda k: expected[k])
    damaged = dict(values)
    damaged[max(values)] += p
    assert agrees_with(damaged, lambda k: expected[k]), 'corruption control failed'
    # Explicit low-order polynomials printed in DLMF 10.41.11.
    low = {0:P(1), 1:(7*p**3-9*p)/QQ(24),
           2:(-135*p**2+594*p**4-455*p**6)/QQ(1152),
           3:(-42525*p**3+451737*p**5-883575*p**7+475475*p**9)/QQ(414720)}
    for k, value in low.items():
        assert values[k] == value
    for k, v in values.items():
        assert v.degree() == 3*k
        assert v(-p) == (-1)**k*v
        assert v.valuation() == k
        if k:
            assert v(0) == 0
    return {'entries':len(values), 'independent_full_polynomial_matches':len(values),
            'dlmf_explicit_matches':len(low), 'controls':'Gaussian moments and deliberate corruption passed',
            'structure':'degree, parity, lowest exponent, origin checked on every applicable entry',
            'measure':measure(values)}

def main():
    if '--stored' in sys.argv:
        path = pathlib.Path(sys.argv[sys.argv.index('--stored')+1])
        document = json.loads(path.read_text())
        block = document['Numbers']
        if isinstance(block, dict):
            values = {int(k):P(v['number'] if isinstance(v,dict) else v) for k,v in block.items()}
        else:
            values = {int(v['params']['k']):P(v['number']) for v in block}
    else:
        from generate import DebyeDerivative
        g = DebyeDerivative()
        params = list(g.enumerate())
        assert params == [{'k':k} for k in range(13)]
        values = {d['k']:P(g.value(d,100)) for d in params}
    assert list(values) == list(range(13))
    assert not names_its_rings(str(pathlib.Path(__file__).with_name('generate.py')))
    report = validate(values)
    if '--companion' in sys.argv:
        document = json.loads(pathlib.Path(sys.argv[sys.argv.index('--companion')+1]).read_text())
        block = document['Numbers']
        if isinstance(block, dict):
            u = {int(k):P(v['number'] if isinstance(v,dict) else v) for k,v in block.items()}
        else:
            u = {int(v['params']['k']):P(v['number']) for v in block}
        assert u[0] == values[0] == 1
        for k in range(1,13):
            assert values[k] == u[k] - p*(1-p**2)*u[k-1]/QQ(2) - p**2*(1-p**2)*u[k-1].derivative()
        assert all(values[k](1) == u[k](1) for k in range(13))
        report['companion_T567'] = 'All 13 endpoint identities and 12 derivative-conversion identities passed'
    print(json.dumps(report,indent=2,default=str))
    print('VALUES_JSON='+json.dumps({str(k):str(v) for k,v in values.items()}))
    return 0

if __name__ == '__main__':
    sys.exit(main())