back to table · edit · history · where entries came from · files · download
9556 bytes, as of the version from 2026-09-27 14:34 (current). Recorded here, not run.
"""100-significant-digit Mathieu values for NumberDB T438.
Install dependencies once and run this standalone file with Python:
python -m pip install numberdb mpmath
python generate.py # verify all existing entries; no writes
python generate.py --sample 10 # verify a smaller sample
python generate.py --publish # explicit value refinement; API key required
Set NUMBERDB_API_KEY and NUMBERDB_ASSISTED_BY in the environment when publishing.
Publishing does not request table promotion, corrections, precision lowering,
entry removal or restatement. This file is self-contained: no repository,
helper files, baseline document, Sage, SciPy or table-value cache is required.
Parameter lists below are copied from the authenticated live snapshot.
Finite Fourier eigensystems use DLMF 28.4 normalization and recurrences:
https://dlmf.nist.gov/28.4 . Sign follows q=0 by continuity (DLMF 28.2).
Both 190-digit/max-harmonic-140 and 240-digit/max-harmonic-180 settings must
support each printed interval before any approximate value is returned.
Each (q, parity) eigensystem pair is cached for all orders/angles in that block.
Rigour is heuristic (agreement-checked), not proven. Finite-matrix residuals
and agreement of larger truncations are not an infinite-operator tail bound.
This generator refuses output precision other than the checked 100-digit target.
"""
import argparse
from decimal import Decimal
from fractions import Fraction
from functools import lru_cache
from pathlib import Path
import mpmath as mp
import numberdb
TABLE = 'T438'
PARAMETERS = ('n', 'q', 'z')
N_VALUES = (1, 2, 3, 4, 5, 6)
Q_VALUES = ('1/8', '1/6', '1/4', '1/3', '1/2', '2/3', '3/4', '1', '4/3', '3/2', '2', '3', '4', '6', '8')
Z_VALUES = ('0', '1/12', '1/10', '1/8', '1/6', '1/5', '1/4', '3/10', '1/3', '3/8', '2/5', '5/12', '1/2')
PROFILES = ((190, 140), (240, 180))
DIAGNOSTICS = {}
def rational(text):
value = Fraction(str(text))
return mp.mpf(value.numerator) / value.denominator
def matrix_for(table, parity, coupling, max_frequency):
first = parity if table == 'T437' else (1 if parity else 2)
harmonics = tuple(range(first, max_frequency + 1, 2))
matrix = mp.matrix(len(harmonics))
for index, harmonic in enumerate(harmonics):
matrix[index, index] = harmonic * harmonic
if harmonic == 1:
matrix[index, index] += coupling if table == 'T437' else -coupling
if index + 1 < len(harmonics):
edge = mp.sqrt(2) * coupling if harmonic == 0 else coupling
matrix[index, index + 1] = edge
matrix[index + 1, index] = edge
return harmonics, matrix
def solution(table, order, harmonics, matrix, eigenvalues, eigenvectors):
rank = order // 2 if table == 'T437' else (order - 1) // 2
coefficients = [eigenvectors[index, rank] for index in range(len(harmonics))]
anchor = mp.fsum(coefficients) if table == 'T437' else mp.fdot(harmonics, coefficients)
if table == 'T437' and harmonics[0] == 0:
anchor += (1 / mp.sqrt(2) - 1) * coefficients[0]
if not mp.isfinite(anchor) or anchor == 0:
raise ArithmeticError('Invalid endpoint sign anchor')
if anchor < 0:
coefficients = [-value for value in coefficients]
eigenvalue = eigenvalues[rank]
norm_error = abs(mp.fdot(coefficients, coefficients) - 1)
residual = max(abs(
mp.fsum(matrix[index, column] * coefficients[column]
for column in range(max(0, index - 1), min(len(harmonics), index + 2)))
- eigenvalue * coefficients[index]) for index in range(len(harmonics)))
tolerance = mp.power(10, -(mp.mp.dps - 20))
if not (mp.isfinite(residual) and mp.isfinite(norm_error)):
raise ArithmeticError('Nonfinite eigenvector diagnostics')
if max(residual, norm_error) > tolerance:
raise ArithmeticError('Finite eigensystem residual or normalization failed')
return coefficients, eigenvalue, residual, norm_error
def evaluate(table, harmonics, coefficients, angle):
argument = mp.pi * rational(angle)
terms = []
for harmonic, coefficient in zip(harmonics, coefficients):
if table == 'T437':
basis = mp.cos(harmonic * argument) if harmonic else 1 / mp.sqrt(2)
else:
basis = mp.sin(harmonic * argument)
terms.append(coefficient * basis)
value = mp.fsum(terms)
if not mp.isfinite(value):
raise ArithmeticError('Nonfinite Fourier sum')
return value
def is_exact_zero(order, angle):
angle = Fraction(angle)
if TABLE == 'T437':
return order % 2 == 1 and angle == Fraction(1, 2)
return angle == 0 or (order % 2 == 0 and angle == Fraction(1, 2))
@lru_cache(maxsize=30)
def checked_group(coupling_text, parity):
if coupling_text not in Q_VALUES or parity not in (0, 1):
raise ValueError('Coupling/parity is outside the fixed live parameter selection')
orders = [order for order in N_VALUES if order % 2 == parity]
results = []
profiles = []
for working_digits, max_frequency in PROFILES:
with mp.workdps(working_digits):
harmonics, matrix = matrix_for(TABLE, parity, rational(coupling_text), max_frequency)
eigenvalues, eigenvectors = mp.eigsy(matrix)
values = {}
diagnostics = {}
for order in orders:
coefficients, eigenvalue, residual, norm_error = solution(
TABLE, order, harmonics, matrix, eigenvalues, eigenvectors)
diagnostics[str(order)] = {'residual': mp.nstr(residual, 25),
'norm_error': mp.nstr(norm_error, 25)}
for angle in Z_VALUES:
values[order, angle] = evaluate(TABLE, harmonics, coefficients, angle)
results.append(values)
profiles.append({'dps': working_digits, 'max_frequency': max_frequency,
'orders': diagnostics})
output = {}
with mp.workdps(PROFILES[-1][0]):
worst_difference = mp.mpf(0)
for order in orders:
for angle in Z_VALUES:
first, second = (values[order, angle] for values in results)
if is_exact_zero(order, angle):
if max(abs(first), abs(second)) > mp.mpf('1e-170'):
raise ArithmeticError('Symmetry zero control failed')
output[order, angle] = 0
continue
text = mp.nstr(first, 100, strip_zeros=False)
if Decimal(text).is_zero():
raise ArithmeticError('Unexpected zero: refuse an approximate exactness claim')
unit = mp.power(10, Decimal(text).as_tuple().exponent)
difference = abs(first - second) / unit
worst_difference = max(worst_difference, difference)
if difference > mp.mpf('0.001'):
raise ArithmeticError('The two settings do not support 100 printed digits')
if abs(mp.mpf(text) - second) > unit:
raise ArithmeticError('Stronger setting is outside the printed interval')
output[order, angle] = text
DIAGNOSTICS[coupling_text, parity] = {
'profiles': profiles, 'worst_profile_difference_in_output_units': mp.nstr(worst_difference, 30)}
return output
class MathieuValues(numberdb.Generator):
table = TABLE
parameters = PARAMETERS
type = 'R'
digits = 100
rigour = 'heuristic (agreement-checked)'
files = (Path(__file__).name,)
def enumerate(self):
for order in N_VALUES:
for coupling in Q_VALUES:
for angle in Z_VALUES:
yield dict(zip(PARAMETERS, (str(order), coupling, angle)))
def value(self, params, digits):
if digits != 100:
raise ValueError('This generator supports the agreement-checked 100-digit target only')
if set(params) != set(PARAMETERS):
raise ValueError('Unexpected parameter names')
order_text, coupling, angle = (str(params[name]) for name in PARAMETERS)
if order_text not in tuple(str(order) for order in N_VALUES):
raise ValueError('Order outside the fixed live parameter selection')
if coupling not in Q_VALUES or angle not in Z_VALUES:
raise ValueError('Coupling/angle outside the fixed live parameter selection')
order = int(order_text)
if is_exact_zero(order, angle):
return 0
return checked_group(coupling, order % 2)[order, angle]
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--publish', action='store_true')
parser.add_argument('--sample', type=int, help='verify this many spread-out entries; default all')
args = parser.parse_args()
if args.sample is not None and args.sample < 1:
parser.error('--sample must be positive')
if args.publish and args.sample is not None:
parser.error('--sample is read-only verification, not partial publication')
generator = MathieuValues()
if args.publish:
print(generator.publish(
message='Refine Mathieu values using two high-precision Fourier settings',
correcting=False, lowering=False, removing=False, restating=False))
return 0
report = generator.verify(sample=args.sample)
print(report)
return 0 if report.ok else 1
if __name__ == '__main__':
raise SystemExit(main())