"""Debye derivative polynomials V_k(p) -- numberdb.org/T568.

Run with SageMath after installing numberdb:
    sage -pip install numberdb
    sage -python generate.py             # verify every stored polynomial
    sage -python checks.py               # independent exact checks, offline
    sage -python generate.py --fill      # fill this draft, API key required

The --fill action writes entries to an existing draft; it does not offer or
publish the table. NUMBERDB_API_KEY authenticates it. The DLMF 10.41.9
normalization retains all rational coefficients in QQ[p], with the signs
outside the K derivative expansion. The default prefix is 0 <= k <= 12.
See argument-selection.md for the finite-range decision and checks.py for
an independent saddle expansion of the differentiated Bessel integral.
"""
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

R = PolynomialRing(QQ, 'p')
p = R.gen()
DEFAULT_MAX_K = 12

@lru_cache(maxsize=None)
def debye_u(k):
    """DLMF differential-integral recurrence with exact coefficient integration."""
    if k == 0:
        return R(1)
    previous = debye_u(k - 1)
    integrand = (1 - 5*p**2)*previous
    integral = R([QQ(0)] + [c/QQ(j + 1) for j, c in enumerate(integrand.list())])
    return p**2*(1-p**2)*previous.derivative()/QQ(2) + integral/QQ(8)

@lru_cache(maxsize=None)
def debye_v(k):
    if k < 0:
        raise ValueError('The correction order must be nonnegative.')
    if k == 0:
        return R(1)
    previous = debye_u(k - 1)
    return (debye_u(k) - p*(1-p**2)*previous/QQ(2)
            - p**2*(1-p**2)*previous.derivative())

class DebyeDerivative(numberdb.Generator):
    table = 'T568'
    parameters = ('k',)
    type = 'Q[]'
    rigour = 'exact'
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
             'numerical-evidence.md', 'skill-manifest.json')

    def enumerate(self, max_k=DEFAULT_MAX_K):
        for k in range(max_k + 1):
            yield {'k': k}

    def value(self, params, digits):
        return debye_v(int(params['k']))

if __name__ == '__main__':
    generator = DebyeDerivative()
    if '--fill' in sys.argv:
        print(generator.publish(message='Fill exact derivative polynomials after independent integral checks'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)