back to table · edit · history · where entries came from · files · download
9117 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Independently check every selected value; never write to the website.
python checks.py T490 --directory .
Checks use finite, narrow Arb enclosures, exact decimal endpoint comparisons,
independent mpmath evaluations at 160 and 210 digits, and elementary integral
representations for the Kelvin and Struve values. The integrals are numerical
agreement checks, not the foundation of the proven error claim.
"""
import argparse
import json
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import mpmath as mp
from flint import acb, arb, ctx
from intrinsic_values import (IntrinsicValues, as_interval, evaluate, mp_derivative, mp_value,
precision, sign, values_and_derivatives)
def written_interval(text):
text = text.rstrip('?')
center = Fraction(text)
unit = Fraction(10) ** Decimal(text).as_tuple().exponent if '.' in text or 'e' in text.lower() else Fraction(0)
return center - unit, center + unit
def contains_ball(text, value):
lower, upper = written_interval(text)
if isinstance(value, int):
return lower <= value <= upper
if not value.is_finite() or value.rel_accuracy_bits() < 350:
return False
enclosed = as_interval(value)
return lower <= enclosed.lower <= enclosed.upper <= upper
def mp_argument(generator, row):
params = row['params']
if 'bracket' in row:
argument, unused = generator.argument(params)
enclosure = as_interval(argument)
midpoint = (enclosure.lower + enclosure.upper) / 2
return mp.mpf(midpoint.numerator) / midpoint.denominator
if row['argument'] == 'turning':
coupling = mp.mpf(params['eta'])
order = int(params['ell'])
return coupling + mp.sqrt(coupling**2 + order * (order + 1))
return mp.mpf(row['argument'])
def independent_integral(kind, order, argument):
if kind in ('ber', 'bei'):
rotated = argument * mp.exp(3j * mp.pi / 4)
result = mp.quad(lambda angle: mp.cos(rotated * mp.sin(angle)),
[0, mp.pi / 4, mp.pi / 2, mp.pi]) / mp.pi
return result.real if kind == 'ber' else result.imag
if kind in ('ker', 'kei'):
scale = argument / mp.sqrt(2)
trig = mp.cos if kind == 'ker' else mp.sin
result = mp.quad(lambda point: mp.exp(-scale * mp.cosh(point)) *
trig(scale * mp.cosh(point)), [0, 1, 2, 4, 8])
return result if kind == 'ker' else -result
trig = mp.sin if kind == 'struveh' else mp.sinh
factor = 2 * (argument / 2)**order / (mp.sqrt(mp.pi) * mp.gamma(order + mp.mpf('0.5')))
return factor * mp.quad(lambda angle: trig(argument * mp.cos(angle)) *
mp.sin(angle)**(2 * order), [0, mp.pi / 4, mp.pi / 2])
def coulomb_identity(params, argument):
order = int(params['ell'])
coupling = arb(params['eta'])
regular = acb(argument).coulomb_f(order, coupling).real
irregular = acb(argument).coulomb_g(order, coupling).real
factor = (1 + coupling**2 / (order + 1)**2).sqrt()
shifted = (order + 1) / argument + coupling / (order + 1)
regular_prime = shifted * regular - factor * acb(argument).coulomb_f(order + 1, coupling).real
irregular_prime = shifted * irregular - factor * acb(argument).coulomb_g(order + 1, coupling).real
wronskian = irregular * regular_prime - regular * irregular_prime
if not wronskian.is_finite() or not wronskian.contains(1) or wronskian.rad() >= arb('1e-110'):
raise ArithmeticError('Failed Coulomb Wronskian')
if coupling == 0:
regular_reference = (arb.pi() * argument / 2).sqrt() * acb(argument).bessel_j(arb(order) + arb('1/2')).real
irregular_reference = -(arb.pi() * argument / 2).sqrt() * acb(argument).bessel_y(arb(order) + arb('1/2')).real
for value, reference in ((regular, regular_reference), (irregular, irregular_reference)):
if not (value - reference).contains(0):
raise ArithmeticError('Failed free-case Bessel identity')
def check(tid, directory, document=None):
generator = IntrinsicValues(tid, directory / 'selection.json')
document = document or json.loads((directory / 'table.json').read_text())
records = []
with precision(900), mp.workdps(150):
if not contains_ball('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068', arb.pi()):
raise ArithmeticError('Decimal interval control failed')
if contains_ball('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117168', arb.pi()):
raise ArithmeticError('Displaced-value negative control failed')
if abs(mp.quad(lambda point: point, [0, 1]) - mp.mpf('0.5')) > mp.mpf('1e-140'):
raise ArithmeticError('Nonzero quadrature control failed')
for row in generator.selection['rows']:
params = row['params']
stored = document['Numbers']
for parameter in generator.parameters:
stored = stored[params[parameter]]
text = stored['number'] if isinstance(stored, dict) else stored
argument, root_proof = generator.argument(params)
value = evaluate(generator.kind, params, argument)
if not contains_ball(text, value):
raise ArithmeticError('Stored interval misses recomputed ball: ' + str(params))
record = {'params': params, 'stored_interval_contains_ball': True,
'accuracy_bits': None if isinstance(value, int) else int(value.rel_accuracy_bits()),
'argument': argument.str(130), 'root_certificate': root_proof}
if root_proof:
record['independent_argument_policy'] = 'Evaluate at the certified bracket midpoint; the rigorous Arb value separately propagates the entire argument ball.'
for decimal_precision in (160, 210):
with mp.workdps(decimal_precision):
independent_argument = mp_argument(generator, row)
if independent_argument == 0:
independent = mp.mpf(1 if generator.kind == 'ber' else 0)
if generator.kind == 'kei':
independent = -mp.pi / 4
elif generator.kind.startswith('coulomb'):
independent = getattr(mp, generator.kind)(int(params['ell']),
mp.mpf(params['eta']), independent_argument)
else:
independent = mp_value(generator.kind, int(params['nu']), independent_argument)
lower, upper = written_interval(text)
if not mp.mpf(lower.numerator) / lower.denominator <= independent <= mp.mpf(upper.numerator) / upper.denominator:
raise ArithmeticError('mpmath value misses stored interval: ' + str(params))
record['mpmath_%d_digits' % decimal_precision] = True
if generator.kind.startswith('coulomb'):
coulomb_identity(params, argument)
if row['argument'] == 'turning':
coefficient = 1 - 2 * arb(params['eta']) / argument - int(params['ell']) * (int(params['ell']) + 1) / argument**2
if not coefficient.is_finite() or not coefficient.contains(0) or coefficient.rad() >= arb('1e-150'):
raise ArithmeticError('Not a certified turning point')
record['wronskian_and_free_case'] = True
elif row['argument'] != '0':
with mp.workdps(150):
quadrature_argument = mp_argument(generator, row)
quadrature = independent_integral(generator.kind, int(params['nu']), quadrature_argument)
reference = mp_value(generator.kind, int(params['nu']), quadrature_argument)
relative = abs((quadrature - reference) / reference)
if not relative < mp.mpf('1e-115'):
raise ArithmeticError('Independent integral failed: ' + str(params))
record['integral_relative_error'] = mp.nstr(relative, 10)
records.append(record)
print(tid, params, 'validated', flush=True)
return {'tid': tid, 'count': len(records), 'all_passed': True,
'controls': ['pi_interval', 'displaced_value_rejected', 'integral_t_equals_half'],
'minimum_accuracy_bits': min(record['accuracy_bits'] for record in records
if record['accuracy_bits'] is not None),
'records': records}
if __name__ == '__main__':
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('tid')
parser.add_argument('--directory', type=Path, default=Path(__file__).parent)
arguments = parser.parse_args()
evidence = check(arguments.tid, arguments.directory)
(arguments.directory / 'checks-evidence.json').write_text(json.dumps(evidence, indent=2) + '\n')