generate.py

back to table · edit · history · where entries came from · files · download

2751 bytes, as of the version from 2026-10-28 01:34 (current). Recorded here, not run.

"""Debye polynomials U_k(p) -- numberdb.org/T567

Run with SageMath:
    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python checks.py                  # independent offline checks
    sage -python generate.py --fill-draft    # NUMBERDB_API_KEY required

The last command writes entries to the existing draft; it does not offer,
review, or publish the table. In the campaign, pipe the key on stdin and set
NUMBERDB_KEY_FROM_STDIN=1. NUMBERDB_PUBLISH=1 selects that same fill action.

DLMF 10.41.9 defines the normalization. The coefficient update is its exact
rational derivative-plus-integral recurrence applied to each monomial.
The entire polynomial is retained, with k=0,...,12 by default. No argument
p is sampled. See argument-selection.md for the extent and conventions.
"""
import os
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()
MAX_K = 12


@lru_cache(maxsize=None)
def debye(k):
    """The kth Debye polynomial, including its exact rational normalization."""
    k = int(k)
    if k < 0:
        raise ValueError('k must be nonnegative')
    if k == 0:
        return R.one()
    previous = debye(k - 1)
    coefficients = [QQ(0)] * (3 * k + 1)
    for j, c in enumerate(previous.list()):
        coefficients[j + 1] += c * QQ((2*j + 1)**2) / QQ(8*(j + 1))
        coefficients[j + 3] -= c * QQ((2*j + 1)*(2*j + 5)) / QQ(8*(j + 3))
    return R(coefficients)


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

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

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


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].strip().isidentifier():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    generator = DebyePolynomials()
    if '--fill-draft' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Add exact Debye polynomials through order twelve'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)