back to table · edit · history · where entries came from · files · download
6474 bytes, as of the version from 2026-10-27 01:01 (current). Recorded here, not run.
"""Reproduce a modest Mathieu reference at q=1, without writing the website.
python3 -m pip install numberdb mpmath
python3 generate.py
The q=1 reference is documented in DLMF 28.3; it is not a uniquely selected
physical coupling. The function arguments 0 and pi/2 are symmetry boundaries.
DLMF 28.2 and 28.4 fix the branch, Fourier recurrence and normalization.
Two finite Fourier sections at 190/240 digits and maximum frequencies
140/180 are compared. This is heuristic agreement, not a rigorous infinite
operator enclosure. Endpoint arguments are evaluated by their exact parity
signs, so no floating-point approximation to pi enters the Fourier evaluation.
The separate checks file uses independent Taylor ODE shooting.
"""
from decimal import Decimal
from functools import lru_cache
import mpmath as mp
import numberdb
PROFILES = ((190, 140), (240, 180))
PARAMETERS = {"T436": ("q", "n"), "T437": ("n", "q", "z_over_pi"),
"T438": ("n", "q", "z")}
ORDERS = {"T436": tuple(range(11)), "T437": tuple(range(6)),
"T438": tuple(range(1, 7))}
ANGLES = ("0", "1/2")
DIAGNOSTICS = {}
def exact_zero(kind, order, angle):
if kind == "ce":
return angle == "1/2" and order % 2 == 1
return angle == "0" or order % 2 == 0
def endpoint(kind, harmonics, coefficients, angle):
terms = []
for harmonic, coefficient in zip(harmonics, coefficients):
if kind == "ce":
if harmonic == 0:
basis = 1 / mp.sqrt(2)
elif angle == "0":
basis = 1
else:
basis = 0 if harmonic % 2 else (-1) ** (harmonic // 2)
elif angle == "0" or harmonic % 2 == 0:
basis = 0
else:
basis = (-1) ** ((harmonic - 1) // 2)
terms.append(coefficient * basis)
return mp.fsum(terms)
@lru_cache(maxsize=None)
def finite_section(kind, parity, digits, maximum):
with mp.workdps(digits):
first = parity if kind == "ce" else (1 if parity else 2)
harmonics = tuple(range(first, maximum + 1, 2))
matrix = mp.matrix(len(harmonics))
for position, harmonic in enumerate(harmonics):
matrix[position, position] = harmonic * harmonic
if harmonic == 1:
matrix[position, position] += 1 if kind == "ce" else -1
if position + 1 < len(harmonics):
edge = mp.sqrt(2) if harmonic == 0 else mp.mpf(1)
matrix[position, position + 1] = edge
matrix[position + 1, position] = edge
eigenvalues, eigenvectors = mp.eigsy(matrix)
output = {}
maximum_order = 10 if kind == "ce" else 6
for order in range(first, maximum_order + 1, 2):
rank = order // 2 if kind == "ce" else (order - 1) // 2
coefficients = [eigenvectors[position, rank] for position in range(len(harmonics))]
anchor = (endpoint(kind, harmonics, coefficients, "0") if kind == "ce"
else mp.fdot(harmonics, coefficients))
if anchor == 0 or not mp.isfinite(anchor):
raise ArithmeticError("Invalid sign anchor")
if anchor < 0:
coefficients = [-value for value in coefficients]
eigenvalue = eigenvalues[rank]
residual = max(abs(mp.fsum(matrix[row, column] * coefficients[column]
for column in range(max(0, row - 1),
min(len(harmonics), row + 2)))
- eigenvalue * coefficients[row])
for row in range(len(harmonics)))
norm_error = abs(mp.fdot(coefficients, coefficients) - 1)
if max(residual, norm_error) > mp.power(10, -digits + 20):
raise ArithmeticError("Finite-section eigenvector check failed")
output[order] = {"eigenvalue": +eigenvalue,
"0": +endpoint(kind, harmonics, coefficients, "0"),
"1/2": +endpoint(kind, harmonics, coefficients, "1/2"),
"residual": +residual, "norm_error": +norm_error}
print("finite section", kind, parity, digits, maximum, "done", flush=True)
return output
def checked_value(tid, params):
if params.get("q") != "1":
raise ValueError("Only the documented q=1 reference is selected")
order = int(params["n"])
if order not in ORDERS[tid]:
raise ValueError("Order outside the selected finite prefix")
kind = "se" if tid == "T438" else "ce"
angle = params.get("z_over_pi", params.get("z"))
key = "eigenvalue" if tid == "T436" else angle
if tid != "T436" and angle not in ANGLES:
raise ValueError("Only symmetry endpoints are selected")
if tid != "T436" and exact_zero(kind, order, angle):
return 0
values = [finite_section(kind, order % 2, *profile)[order][key] for profile in PROFILES]
with mp.workdps(PROFILES[-1][0]):
if not all(mp.isfinite(value) for value in values):
raise ArithmeticError("Nonfinite Fourier value")
written = mp.nstr(values[-1], 100, strip_zeros=False)
unit = mp.power(10, Decimal(written).as_tuple().exponent)
discrepancy = abs(values[0] - values[1]) / unit
if discrepancy >= mp.mpf("0.001"):
raise ArithmeticError("Finite sections disagree at the printed precision")
DIAGNOSTICS[tid + ":" + ",".join(params[name] for name in PARAMETERS[tid])] = {
"difference_ulps": mp.nstr(discrepancy, 20),
"profiles": PROFILES,
}
return written
class MathieuReference(numberdb.Generator):
type = "R"
digits = 100
rigour = "heuristic (agreement-checked)"
def __init__(self, table):
if table not in PARAMETERS:
raise ValueError("Unassigned Mathieu table")
self.table = table
self.parameters = PARAMETERS[table]
def enumerate(self):
for order in ORDERS[self.table]:
if self.table == "T436":
yield {"q": "1", "n": str(order)}
else:
for angle in ANGLES:
yield dict(zip(self.parameters, (str(order), "1", angle)))
def value(self, params, digits):
if digits != 100:
raise ValueError("Only the validated 100-digit target is supported")
return checked_value(self.table, params)