back to table · edit · history · where entries came from · files · download
5559 bytes, as of the version from 2026-10-28 01:47 (current). Recorded here, not run.
"""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())