back to table · edit · history · where entries came from · files · download
5998 bytes, as of the version from 2026-10-28 01:34 (current). Recorded here, not run.
"""Independent exact checks for Debye polynomials -- numberdb.org/T567.
sage -python checks.py
sage -python checks.py --stored # authenticate for an unpublished draft
The independent computation expands the Bessel K integral by Gaussian
moments. It does not use the Debye recurrence. The derivation is attached as
independent-derivation.md. The checks also use reciprocal-Stirling
coefficients from Bernoulli numbers and three published polynomial controls.
"""
import json
import os
import sys
from math import factorial, comb
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 check import exactness, measure, agrees_with, names_its_rings
P = PolynomialRing(QQ, 'p')
p = P.gen()
A = PolynomialRing(QQ, names=('a', 'v'))
a, v = A.gens()
def gaussian_moment(n):
if n % 2:
return 0
answer = 1
for j in range(1, n, 2):
answer *= j
return answer
def saddle_polynomials(max_k):
# a=sqrt(p), u=a*v/sqrt(nu); h=1/sqrt(nu).
assert [gaussian_moment(j) for j in range(7)] == [1, 0, 1, 0, 3, 0, 15]
exponent = [A.zero()]
for j in range(1, 2*max_k+1):
power = j if j % 2 == 0 else j+2
exponent.append(-a**power * v**(j+2) / QQ(factorial(j+2)))
expanded = [A.one()]
answer = []
for n in range(2*max_k+1):
if n:
expanded.append(sum((j*exponent[j]*expanded[n-j]
for j in range(1, n+1)), A.zero()) / QQ(n))
integrated = P.zero()
for (a_power, v_power), c in expanded[n].dict().items():
moment = gaussian_moment(v_power)
if moment:
assert a_power % 2 == 0
integrated += c*moment*p**(a_power//2)
if n % 2:
assert integrated == 0
else:
answer.append((-1)**(n//2)*integrated)
assert answer[0] == 1
assert answer[1] == (3*p-5*p**3)/QQ(24)
return answer
def reciprocal_stirling(max_k):
bernoulli = [QQ(1)]
for n in range(1, max_k+2):
bernoulli.append(-sum((QQ(comb(n+1,j))*bernoulli[j]
for j in range(n)), QQ(0))/QQ(n+1))
assert bernoulli[1:5] == [QQ(-1)/2, QQ(1)/6, QQ(0), -QQ(1)/30]
log_coefficients = [QQ(0)]*(max_k+1)
for j in range(1,max_k+1,2):
log_coefficients[j] = -bernoulli[j+1]/QQ(j*(j+1))
coefficients = [QQ(1)]
for n in range(1,max_k+1):
coefficients.append(sum((j*log_coefficients[j]*coefficients[n-j]
for j in range(1,n+1)), QQ(0))/QQ(n))
assert coefficients[:4] == [QQ(1),-QQ(1)/12,QQ(1)/288,QQ(139)/51840]
return coefficients
def validate(values):
assert list(values) == list(range(13)), 'Expected every k=0,...,12 in order'
assert not exactness(values)
independent = saddle_polynomials(12)
assert not agrees_with(values, lambda k: independent[k])
controls = [P(1), (3*p-5*p**3)/QQ(24),
(81*p**2-462*p**4+385*p**6)/QQ(1152),
(30375*p**3-369603*p**5+765765*p**7-425425*p**9)/QQ(414720)]
assert [values[k] for k in range(4)] == controls
gamma = reciprocal_stirling(12)
for k, u in values.items():
assert u.degree() == 3*k
assert u(-p) == (-1)**k*u
assert u.valuation() == k
assert u(1) == gamma[k]
assert u(0) == (1 if k == 0 else 0)
if k < 12:
integrand = (1-5*p**2)*u
integral = sum((c*p**(j+1)/QQ(j+1)
for j,c in enumerate(integrand.list())), P.zero())
assert values[k+1] == p**2*(1-p**2)*u.derivative()/QQ(2)+integral/QQ(8)
# A deliberately wrong constant and coefficient must fail the independent comparison.
corrupted = dict(values); corrupted[12] += p**12
assert agrees_with(corrupted, lambda k: independent[k])
return {'exactness':'passed', 'saddle_full_polynomials':13,
'published_controls_k':'0..3', 'reciprocal_stirling_endpoints':13,
'recurrence_transitions':12, 'degree_parity_valuation_origin':13,
'gaussian_moment_control':[1,0,1,0,3,0,15],
'deliberate_corruption_detected':True, 'measure':measure(values),
'values':{str(k):str(u) for k,u in values.items()}}
def stored_values(document):
"""Read this one-parameter table in either API storage representation."""
entries = document['Numbers']
if isinstance(entries, list):
return {int(e['params']['k']): P(e['number']) for e in entries}
assert isinstance(entries, dict), 'Expected flat or nested number records'
return {int(k): P(e['number'] if isinstance(e, dict) else e)
for k, e in entries.items()}
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].strip().isidentifier():
key=key.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY']=key
stored = '--stored' in sys.argv or os.environ.get('NUMBERDB_RESTATING') == '1'
if stored:
document=numberdb.table('T567')
# The API returns the document's number strings. Parse them into QQ[p].
values=stored_values(document)
mode='stored'
else:
from generate import DebyePolynomials
generator=DebyePolynomials()
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
values={int(params['k']):generator.value(params,100)
for params in generator.enumerate()}
mode='offline'
evidence=validate(values)
evidence['mode']=mode
from sage.version import version as sage_version
evidence['versions']={'sage':sage_version,'python':sys.version.split()[0],
'numberdb':getattr(numberdb,'__version__','unknown')}
print(json.dumps(evidence,indent=2))
if __name__ == '__main__':
main()