"""Independent exact checks for Hankel coefficients -- numberdb.org/T569.

Run `sage -python checks.py` for an offline check, or set
NUMBERDB_CHECK_STORED=1 to check the API readback and the stored Bessel
polynomials. NUMBERDB_KEY_FROM_STDIN=1 authenticates private-draft reads.

The independent reconstruction uses y_0=1, y_1=x+1 and
 y_n=(2*n-1)*x*y_(n-1)+y_(n-2), then interpolates the coefficient of x^k
at t=(n+1/2)^2. DLMF 10.49.1 and 10.49.12 identify this coefficient with
a_k(n+1/2). Evenness and degree 2k turn 21 values into a polynomial
identity for each k<=20. No product or k-recurrence is used to build the
interpolant. The product recurrence and modified-Bessel ODE are additional
consistency checks, not claimed as independent computations.
"""
import os
import sys
import json
from math import factorial
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.env import SAGE_VERSION
from generate import HankelCoefficients, MAX_K, R, nu
from check import exactness, measure, agrees_with, names_its_rings

P = PolynomialRing(QQ, 'x')
x = P.gen()
S = PolynomialRing(QQ, 't')
t = S.gen()


def bessel_polynomials(max_n):
    b = [P.one(), 1+x]
    for n in range(2, max_n+1):
        b.append((2*n-1)*x*b[-1]+b[-2])
    assert b[0] == 1 and b[1] == 1+x
    assert b[2] == 1+3*x+3*x**2
    return b[:max_n+1]


def interpolate_from_bessel(b, max_k):
    nodes = [(QQ(n)+QQ(1)/2)**2 for n in range(max_k+1)]
    basis = []
    for n, node in enumerate(nodes):
        term = S.one()
        for m, other in enumerate(nodes):
            if m != n:
                term *= (t-other) * (QQ(1)/(node-other))
        basis.append(term)
    # Nonzero controls: interpolation preserves 1, t and t^2 exactly.
    for exponent in range(3):
        assert sum(nodes[n]**exponent*basis[n] for n in range(max_k+1)) == t**exponent
    return {k: R(sum(b[n][k]*basis[n] for n in range(max_k+1))(nu**2))
            for k in range(max_k+1)}


def verify_values(values, check_live_bessel=False):
    assert list(values) == list(range(MAX_K+1)), 'unexpected enumeration/order'
    assert not exactness(values)
    assert not names_its_rings(Path(__file__).with_name('generate.py'))
    b = bessel_polynomials(MAX_K)
    independent = interpolate_from_bessel(b, MAX_K)
    assert not agrees_with(values, lambda k: independent[k])
    corrupt = dict(values)
    corrupt[MAX_K] += QQ(1)/7
    assert len(agrees_with(corrupt, lambda k: independent[k])) == 1
    specializations = roots = 0
    for k,a in values.items():
        assert a.degree() == 2*k and a(-nu) == a
        assert a.leading_coefficient() == QQ(1)/QQ(2**k*factorial(k))
        if k:
            assert QQ(8*k)*a == (4*nu**2-(2*k-1)**2)*values[k-1]
        else:
            assert a == 1
        for n in range(MAX_K+1):
            expected = (QQ(factorial(n+k))/QQ(2**k*factorial(k)*factorial(n-k))
                        if k<=n else QQ(0))
            assert a(QQ(n)+QQ(1)/2) == b[n][k] == expected
            specializations += 1
        for j in range(1,k+1):
            assert a(QQ(2*j-1)/2) == a(-QQ(2*j-1)/2) == 0
            roots += 2
    # t=1/z; remove the exponential and square-root prefactor from K_nu.
    T = PolynomialRing(R, 't')
    u = T.gen()
    f = sum(T(a)*u**k for k,a in values.items())
    residual = u**2*f.derivative(2)+(2*u+2)*f.derivative()+(QQ(1)/4-nu**2)*f
    assert all(residual[j] == 0 for j in range(MAX_K))
    assert residual[MAX_K] == (QQ(MAX_K*(MAX_K+1))+QQ(1)/4-nu**2)*values[MAX_K]
    live_count = 0
    if check_live_bessel:
        doc = numberdb.table('Bessel_polynomials')
        for n in range(MAX_K+1):
            raw = doc['Numbers'][str(n)]
            if isinstance(raw,dict): raw=raw['number']
            assert P(raw) == b[n]
            live_count += 1
    return {
        'all_passed': True, 'exact_polynomials': len(values),
        'independent_method': 'exact interpolation in nu^2 from the Bessel polynomial n-recurrence',
        'independent_polynomial_equalities': len(values),
        'half_integer_comparisons': specializations,
        'roots_checked_with_both_signs': roots,
        'recurrence_checks': MAX_K,
        'ode_coefficients_vanishing': MAX_K,
        'stored_bessel_polynomials_compared': live_count,
        'controls': ['y_0=1, y_1=x+1, y_2=3*x^2+3*x+1',
                     'interpolation reproduces 1,t,t^2',
                     'adding 1/7 to a_20 is detected'],
        'measurements': measure(values),
        'sage_version': SAGE_VERSION,
        'skill_manifest': 'agents/runs/20261028T022135Z-build-skills/manifest.json',
        'skill_sha256': 'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c'
    }


def get_stored_values():
    doc=numberdb.table('T569')
    block=doc['Numbers']
    if isinstance(block,dict):
        return {int(k):R(v['number'] if isinstance(v,dict) else v) for k,v in block.items()}
    return {int(v['params']['k']):R(v['number']) for v in block}


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key=sys.stdin.read().strip()
        if '=' in key and key.split('=',1)[0].replace('_','').isalnum():
            key=key.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY']=key
    stored=os.environ.get('NUMBERDB_CHECK_STORED') == '1'
    g=HankelCoefficients()
    values=(get_stored_values() if stored else
            {int(p['k']):g.value(p,100) for p in g.enumerate()})
    report=verify_values(values, check_live_bessel=True)
    report['values_read_back_from_api']=stored
    report['values']={str(k):str(v) for k,v in values.items()}
    print(json.dumps(report,indent=2))