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()
def debye(k):
U = R(1)
for _ in range(k):
f = (1 - 5*p**2)*U
integral = R([QQ(0)] + [c/QQ(j+1) for j, c in enumerate(f.list())])
U = p**2*(1-p**2)*U.derivative()/QQ(2) + integral/QQ(8)
return U
values = [debye(k) for k in range(13)]
values
All coefficients were computed exactly in $\mathbb{Q}[p]$ from (1). An independent expansion of the Bessel integral [3], using exact Gaussian moments, reproduces every coefficient of all thirteen polynomials. The first three nonconstant polynomials agree with DLMF 10.41.10.
The reciprocal scaled-Gamma expansion [2] independently reproduces all thirteen values at $p=1$. The recurrence, degree, parity and lowest exponent were checked over the stored range. The attached checks include known Gaussian moments and a deliberate coefficient corruption as controls.