back to table · edit · history · where entries came from · files · download
9927 bytes, as of the version from 2026-09-20 18:00 (current). Recorded here, not run.
"""Secondary polynomials of the Hermite polynomials q_n -- numberdb.org/T373.
For the physicists' Hermite polynomial H_n, this table stores
q_n(x) = pi^(-1/2) integral_{-infinity}^{infinity}
(H_n(t) - H_n(x)) / (t - x) exp(-t^2) dt.
Run it with SageMath:
$ sage -pip install numberdb # once
$ sage -python generate.py # check the table against this code
$ sage -python generate.py --publish # fill the draft, with NUMBERDB_API_KEY set
The values are exact rational polynomials. The integrity checks compare the
recurrence with the defining moment formula, with a plain Python fractions
implementation of the same formula, with the shifted monic associated-Hermite
recurrence, and with the Gauss-Hermite quadrature identity
q_n(x_k)/H_n'(x_k) = w_k/sqrt(pi) at the roots of H_n, checked exactly as a
polynomial congruence modulo H_n. They also check the Pade property of
q_n/H_n for the Stieltjes transform of the probability Hermite weight.
"""
import os
import sys
from fractions import Fraction
import numberdb.sage as numberdb
from sage.arith.misc import factorial
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ
TABLE = os.environ.get("NUMBERDB_TABLE", "T373")
UP_TO = 50
R = PolynomialRing(QQ, "x")
x = R.gen()
def _key_from_stdin():
if os.environ.get("NUMBERDB_KEY_FROM_STDIN") != "1":
return
token = sys.stdin.read().strip()
if "=" in token and token.split("=", 1)[0].isupper():
token = token.split("=", 1)[1].strip().strip("'\"")
if token:
os.environ["NUMBERDB_API_KEY"] = token
numberdb.configure(api_key=token)
def hermite_values(up_to=UP_TO):
values = [R.one()]
if up_to == 0:
return values
values.append(ZZ(2) * x)
for n in range(1, up_to):
values.append(ZZ(2) * x * values[n] - ZZ(2 * n) * values[n - 1])
return values
HERMITE = hermite_values(UP_TO)
def secondary_values(up_to=UP_TO):
values = [R.zero()]
if up_to == 0:
return values
values.append(R(2))
for n in range(1, up_to):
values.append(ZZ(2) * x * values[n] - ZZ(2 * n) * values[n - 1])
return values
SECONDARY = secondary_values(UP_TO)
def moment(power):
power = int(power)
if power % 2:
return QQ(0)
k = power // 2
return QQ(factorial(2 * k)) / (QQ(4) ** k * QQ(factorial(k)))
def secondary_from_moments(n):
total = R.zero()
polynomial = HERMITE[int(n)]
for degree, coefficient in enumerate(polynomial.list()):
for i in range(degree):
total += coefficient * moment(i) * x ** (degree - 1 - i)
return R(total)
def fraction_moment(power):
power = int(power)
if power % 2:
return Fraction(0)
k = power // 2
return Fraction(int(factorial(2 * k)), (4 ** k) * int(factorial(k)))
def fraction_poly_add(one, other):
length = max(len(one), len(other))
out = [Fraction(0) for _ in range(length)]
for i, coefficient in enumerate(one):
out[i] += coefficient
for i, coefficient in enumerate(other):
out[i] += coefficient
return out
def fraction_poly_mul_x(poly):
return [Fraction(0)] + list(poly)
def fraction_poly_scalar(poly, scalar):
return [scalar * coefficient for coefficient in poly]
def fraction_hermite_values(up_to=UP_TO):
values = [[Fraction(1)]]
if up_to == 0:
return values
values.append([Fraction(0), Fraction(2)])
for n in range(1, up_to):
values.append(fraction_poly_add(
fraction_poly_scalar(fraction_poly_mul_x(values[n]), Fraction(2)),
fraction_poly_scalar(values[n - 1], Fraction(-2 * n)),
))
return values
FRACTION_HERMITE = fraction_hermite_values(UP_TO)
def fraction_secondary_from_moments(n):
out = [Fraction(0) for _ in range(int(n))]
for degree, coefficient in enumerate(FRACTION_HERMITE[int(n)]):
for i in range(degree):
out[degree - 1 - i] += coefficient * fraction_moment(i)
while len(out) > 1 and out[-1] == 0:
out.pop()
return out
def sage_poly_from_fractions(coefficients):
total = R.zero()
for degree, coefficient in enumerate(coefficients):
total += QQ(coefficient.numerator) / QQ(coefficient.denominator) * x ** degree
return total
def check_identities():
for n in range(UP_TO + 1):
value = SECONDARY[n]
if value != secondary_from_moments(n):
raise ArithmeticError("moment formula failed at n=%d" % n)
if value != sage_poly_from_fractions(fraction_secondary_from_moments(n)):
raise ArithmeticError("plain Python fraction check failed at n=%d" % n)
if n == 0:
if value != 0:
raise ArithmeticError("q_0 is not zero")
continue
if value.degree() != n - 1:
raise ArithmeticError("degree failed at n=%d" % n)
if value.leading_coefficient() != QQ(2) ** n:
raise ArithmeticError("leading coefficient failed at n=%d" % n)
if value(-x) != (-1) ** (n - 1) * value:
raise ArithmeticError("parity failed at n=%d" % n)
for n in range(0, UP_TO):
associated = SECONDARY[n + 1] / (QQ(2) ** (n + 1))
if associated.leading_coefficient() != 1:
raise ArithmeticError("associated polynomial is not monic at n=%d" % n)
if n == 0 and associated != 1:
raise ArithmeticError("associated initial value failed")
if n == 1 and associated != x:
raise ArithmeticError("associated first-degree value failed")
if n >= 1:
previous = SECONDARY[n] / (QQ(2) ** n)
before_previous = SECONDARY[n - 1] / (QQ(2) ** (n - 1))
expected = x * previous - QQ(n) / QQ(2) * before_previous
if associated != expected:
raise ArithmeticError("associated recurrence failed at n=%d" % n)
for n in range(1, UP_TO + 1):
congruence = (SECONDARY[n] * HERMITE[n].derivative()) % HERMITE[n]
expected = R(QQ(2) ** (n + 1) * QQ(factorial(n)))
if congruence != expected:
raise ArithmeticError("Gauss-Hermite congruence failed at n=%d" % n)
for n in range(1, UP_TO + 1):
coefficients = HERMITE[n].list()
for exponent in range(n):
coefficient = sum(
coefficients[degree] * moment(degree + exponent)
for degree in range(len(coefficients))
)
if coefficient != 0:
raise ArithmeticError(
"Pade remainder coefficient %d failed at n=%d"
% (exponent, n)
)
first_remainder = sum(
coefficients[degree] * moment(degree + n)
for degree in range(len(coefficients))
)
if first_remainder == 0:
raise ArithmeticError("first Pade remainder vanished at n=%d" % n)
longest = max((len(str(SECONDARY[n])), n) for n in range(UP_TO + 1))
print("integrity checks passed for q_0 through q_%d" % UP_TO)
print("moment formula and plain Python fraction checks passed")
print("associated-Hermite recurrence and Gauss-Hermite congruence passed")
print("Stieltjes-transform Pade check passed")
print("longest polynomial has %d characters at n=%d" % (longest[0], longest[1]))
class SecondaryHermitePolynomials(numberdb.Generator):
"""Generator for T373, secondary polynomials of physicists' Hermite H_n."""
table = TABLE
parameters = ("n",)
type = "Q[]"
rigour = "exact"
def enumerate(self, up_to=UP_TO):
for n in range(up_to + 1):
yield {"n": str(n)}
def value(self, params, digits):
return SECONDARY[int(params["n"])]
def fill_draft_once(generator, message):
"""Fill a fresh draft without the client's empty upsert probe."""
from numberdb._generate import (
_check_precision,
_check_rigour,
_producer,
_run_name,
_source_files,
)
from numberdb._write import Entries, attach, submit_entries, to_text
table = generator.table
run = _run_name(generator)
entries = Entries(*generator.parameters)
for params in generator.enumerate():
params = dict(params)
wanted = generator.digits_for(params)
entry = generator._entry(params, wanted)
value = entry["number"]
identity = ",".join(str(params[name]) for name in generator.parameters)
_check_rigour(generator, table, identity, value)
written = to_text(value, wanted, generator.format)
_check_precision(table, identity, written, wanted, lowering=False)
record = dict(entry)
record.pop("digits", None)
entries.add(**params, **record, digits=wanted)
answer = submit_entries(
table,
entries,
message=message,
produced_by=_producer(generator, os.environ.get("NUMBERDB_ASSISTED_BY", "")),
upsert=False,
run=run,
rigour=generator.rigour,
)
for name, body in sorted(_source_files(generator).items()):
attach(table, name, body, run=run, message=message,
rigour=generator.rigour)
return answer
if __name__ == "__main__":
_key_from_stdin()
check_identities()
generator = SecondaryHermitePolynomials()
if "--publish" in sys.argv or bool(int(os.environ.get("NUMBERDB_PUBLISH", "0"))):
print(fill_draft_once(
generator,
message="exact Hermite secondary polynomials from recurrence"))
elif os.environ.get("NUMBERDB_API_KEY"):
report = generator.verify(sample=None)
print(report)
if not report.ok:
sys.exit(1)
else:
print("identity checks passed; NUMBERDB_API_KEY is not set, so verify() was skipped")