back to table · edit · history · where entries came from · files · download
14672 bytes, as of the version from 2026-10-27 01:01 (current). Recorded here, not run.
"""Reproduce bounded Mathieu references, without writing the website.
python -m pip install numberdb mpmath
python generate.py --table T436
python generate.py --all --output checks.json
The second command recomputes this table; the last independently checks all
three tables by Taylor ODE shooting. No command publishes or approves data.
The default table is replaced when this file is attached to another table.
The reference couplings 1/2, 1, 2 span modest Fourier coupling strengths in
w'' + (a - 2 q cos(2z)) w = 0. They are a bounded exploratory selection,
not three uniquely preferred physical constants. Orders are the natural
spectral indices, bounded at 10 for characteristic values and at 5/6 for
ce/se. The boundaries 0 and pi/2 and the potential's zero pi/4 provide an
interior reference without an angular mesh. Existing q=1 exact zeros are
retained, but new identically zero boundary values are not enumerated.
DLMF 28.2 and 28.4 define the equation, normalization and Fourier recurrence;
28.6 describes small-q expansions. These sources justify the structure, not
our finite exploratory coupling selection. The Fourier calculations at
160/210 digits and maximum frequencies 110/150 and the independently
implemented Taylor checks are numerical agreement, not infinite-dimensional
operator or ODE remainder enclosures. Arguments are recomputed at each
precision. Endpoint trigonometric factors are exact parity signs.
"""
import argparse
import json
from decimal import Decimal
from functools import lru_cache
from pathlib import Path
import mpmath as mp
import numberdb
TABLE = "T436"
PROFILES = ((160, 110), (210, 150))
ODE_PROFILES = ((160, 12, 140), (190, 16, 170))
COUPLINGS = ("1/2", "1", "2")
PARAMETERS = {"T436": ("q", "n"), "T437": ("n", "q", "z_over_pi"),
"T438": ("n", "q", "z")}
ORDERS = {"T436": range(11), "T437": range(6), "T438": range(1, 7)}
ANGLES = ("0", "1/4", "1/2")
def rational(text):
numerator, separator, denominator = str(text).partition("/")
return mp.mpf(numerator) / (mp.mpf(denominator) if separator else 1)
def exact_zero(kind, order, angle):
return ((kind == "ce" and angle == "1/2" and order % 2 == 1)
or (kind == "se" and (angle == "0"
or (angle == "1/2" and order % 2 == 0))))
def evaluation(kind, harmonics, coefficients, angle):
terms = []
for harmonic, coefficient in zip(harmonics, coefficients):
if kind == "ce" and harmonic == 0:
basis = 1 / mp.sqrt(2)
elif angle == "0":
basis = 1 if kind == "ce" else 0
elif angle == "1/2":
if kind == "ce":
basis = 0 if harmonic % 2 else (-1) ** (harmonic // 2)
else:
basis = (-1) ** ((harmonic - 1) // 2) if harmonic % 2 else 0
else:
phase = harmonic * rational(angle)
basis = mp.cospi(phase) if kind == "ce" else mp.sinpi(phase)
terms.append(coefficient * basis)
return mp.fsum(terms)
@lru_cache(maxsize=None)
def finite_section(kind, parity, coupling_text, digits, maximum):
with mp.workdps(digits):
coupling = rational(coupling_text)
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] += coupling if kind == "ce" else -coupling
if position + 1 < len(harmonics):
edge = coupling * (mp.sqrt(2) if harmonic == 0 else 1)
matrix[position, position + 1] = edge
matrix[position + 1, position] = edge
eigenvalues, eigenvectors = mp.eigsy(matrix)
output = {}
for order in range(first, 11 if kind == "ce" else 7, 2):
rank = order // 2 if kind == "ce" else (order - 1) // 2
coefficients = [eigenvectors[position, rank] for position in range(len(harmonics))]
anchor = (evaluation(kind, harmonics, coefficients, "0") if kind == "ce"
else mp.fdot(harmonics, coefficients))
if not mp.isfinite(anchor) or anchor == 0:
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 normalization or residual failed")
output[order] = {"eigenvalue": +eigenvalue,
**{angle: +evaluation(kind, harmonics, coefficients, angle)
for angle in ANGLES},
"residual": +residual, "norm_error": +norm_error,
"last_coefficient": +abs(coefficients[-1])}
print("Fourier", kind, parity, coupling_text, digits, maximum, flush=True)
return output
def checked_value(tid, params):
order = int(params["n"])
if params["q"] not in COUPLINGS or order not in ORDERS[tid]:
raise ValueError("Outside the validated reference selection")
kind = "se" if tid == "T438" else "ce"
key = "eigenvalue" if tid == "T436" else params[PARAMETERS[tid][-1]]
if tid != "T436" and key not in ANGLES:
raise ValueError("Outside the validated angular selection")
if tid != "T436" and exact_zero(kind, order, key):
return 0
values = [finite_section(kind, order % 2, params["q"], *profile)[order][key]
for profile in PROFILES]
with mp.workdps(PROFILES[-1][0]):
written = mp.nstr(values[-1], 100, strip_zeros=False)
unit = mp.power(10, Decimal(written).as_tuple().exponent)
if not all(mp.isfinite(value) for value in values):
raise ArithmeticError("Nonfinite Fourier value")
if abs(values[0] - values[1]) / unit >= mp.mpf("0.001"):
raise ArithmeticError("Fourier sections disagree at the written precision")
return written
class MathieuReference(numberdb.Generator):
type = "R"
digits = 100
rigour = "heuristic (agreement-checked)"
def __init__(self, table=TABLE):
self.table = table
self.parameters = PARAMETERS[table]
def enumerate(self):
for coupling in COUPLINGS:
for order in ORDERS[self.table]:
if self.table == "T436":
yield {"q": coupling, "n": str(order)}
else:
for angle in ANGLES:
kind = "se" if self.table == "T438" else "ce"
if coupling != "1" and exact_zero(kind, order, angle):
continue
yield dict(zip(self.parameters, (str(order), coupling, 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)
def shoot(kind, eigenvalue, coupling, digits, steps, degree):
with mp.workdps(digits):
step = mp.pi / (2 * steps)
position = mp.mpf(0)
state = [mp.mpf(kind == "ce"), mp.mpf(kind == "se"), mp.mpf(0), mp.mpf(0)]
midpoint = None
for interval in range(steps):
cosine = [mp.cos(2 * position), -2 * mp.sin(2 * position)]
for index in range(2, degree):
cosine.append(-4 * cosine[index - 2] / (index * (index - 1)))
potential = [2 * coupling * value for value in cosine]
potential[0] -= eigenvalue
solution = state[:2]
variation = state[2:]
for index in range(degree - 1):
denominator = (index + 1) * (index + 2)
solution.append(mp.fsum(potential[column] * solution[index - column]
for column in range(index + 1)) / denominator)
variation.append((mp.fsum(potential[column] * variation[index - column]
for column in range(index + 1))
- solution[index]) / denominator)
state = [mp.polyval(list(reversed(solution)), step),
mp.polyval(list(reversed([index * solution[index]
for index in range(1, len(solution))])), step),
mp.polyval(list(reversed(variation)), step),
mp.polyval(list(reversed([index * variation[index]
for index in range(1, len(variation))])), step)]
position += step
if interval + 1 == steps // 2:
midpoint = state[0]
return state, midpoint
def normalized(kind, order, coupling_text, initial, profile):
digits, steps, degree = profile
with mp.workdps(digits):
neumann = order % 2 == (0 if kind == "ce" else 1)
eigenvalue = +initial
coupling = rational(coupling_text)
for iteration in range(2):
state, midpoint = shoot(kind, eigenvalue, coupling, digits, steps, degree)
value, slope, variation, variation_slope = state
eigenvalue -= (slope / variation_slope if neumann else value / variation)
state, midpoint = shoot(kind, eigenvalue, coupling, digits, steps, degree)
value, slope, variation, variation_slope = state
norm = slope * variation - value * variation_slope
if not mp.isfinite(norm) or norm <= 0:
raise ArithmeticError("Invalid independent norm")
scale = mp.sqrt(mp.pi / (4 * norm))
return {"eigenvalue": +eigenvalue, "0": +scale if kind == "ce" else mp.mpf(0),
"1/4": +(scale * midpoint),
"1/2": +(scale * value) if neumann else mp.mpf(0),
"boundary_residual": +(slope if neumann else value)}
def controls():
with mp.workdps(170):
for kind in ("ce", "se"):
for order in (1, 2, 3):
state, midpoint = shoot(kind, mp.mpf(order * order), mp.mpf(0), 160, 12, 140)
value, slope, variation, variation_slope = state
expected_value = mp.cospi(mp.mpf(order) / 2) if kind == "ce" else mp.sinpi(mp.mpf(order) / 2) / order
expected_slope = -order * mp.sinpi(mp.mpf(order) / 2) if kind == "ce" else mp.cospi(mp.mpf(order) / 2)
expected_midpoint = mp.cospi(mp.mpf(order) / 4) if kind == "ce" else mp.sinpi(mp.mpf(order) / 4) / order
expected_norm = mp.pi / 4 if kind == "ce" else mp.pi / (4 * order * order)
if max(abs(value - expected_value), abs(slope - expected_slope),
abs(midpoint - expected_midpoint),
abs(slope * variation - value * variation_slope - expected_norm)) > mp.mpf("1e-130"):
raise ArithmeticError("Zero-coupling ODE control failed")
print("Exact zero-coupling controls passed", flush=True)
def validate(output):
controls()
independent = {}
result = {"fourier_profiles": PROFILES, "ode_profiles": ODE_PROFILES,
"rigour": "heuristic (agreement-checked)", "controls_passed": True,
"tables": {}}
for coupling in COUPLINGS:
for kind in ("ce", "se"):
for order in (range(11) if kind == "ce" else range(1, 7)):
initial = finite_section(kind, order % 2, coupling, *PROFILES[-1])[order]["eigenvalue"]
profiles = [normalized(kind, order, coupling, initial, profile) for profile in ODE_PROFILES]
with mp.workdps(210):
if any(abs(profiles[0][key] - profiles[1][key]) > mp.mpf("1e-125")
for key in ("eigenvalue", *ANGLES)):
raise ArithmeticError("Independent ODE profiles disagree")
if abs(profiles[-1]["boundary_residual"]) > mp.mpf("1e-125"):
raise ArithmeticError("Independent boundary residual too large")
independent[kind, order, coupling] = profiles[-1]
print("ODE", kind, order, coupling, flush=True)
for tid in PARAMETERS:
rows = []
generator = MathieuReference(tid)
for params in generator.enumerate():
written = str(generator.value(params, 100))
kind = "se" if tid == "T438" else "ce"
key = "eigenvalue" if tid == "T436" else params[PARAMETERS[tid][-1]]
with mp.workdps(210):
expected = independent[kind, int(params["n"]), params["q"]][key]
if written == "0":
if expected != 0:
raise ArithmeticError("Exact-zero mismatch")
discrepancy = mp.mpf(0)
else:
unit = mp.power(10, Decimal(written).as_tuple().exponent)
discrepancy = abs(expected - mp.mpf(written)) / unit
if discrepancy >= 1:
raise ArithmeticError("Independent value outside written interval")
rows.append({"parameters": params, "number": written,
"independent_error_ulps": mp.nstr(discrepancy, 20)})
result["tables"][tid] = rows
output.parent.mkdir(parents=True, exist_ok=True)
output.write_text(json.dumps(result, indent=2) + "\n")
print("Independent checks saved", output, flush=True)
if __name__ == "__main__":
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--table", choices=PARAMETERS, default=TABLE)
parser.add_argument("--all", action="store_true")
parser.add_argument("--output", type=Path, default=Path("checks.json"))
arguments = parser.parse_args()
if arguments.all:
validate(arguments.output)
else:
generator = MathieuReference(arguments.table)
for parameters in generator.enumerate():
print(parameters, generator.value(parameters, 100), flush=True)