back to table · edit · history · where entries came from · files · download
4975 bytes, as of the version from 2026-10-27 01:14 (current). Recorded here, not run.
"""Recheck T493's unchanged numbers during its presentation repair.
Run with SageMath and a JSON table document obtained from GET /api/table?id=T493:
sage -pip install numberdb
sage -python repair_checks.py < T493-document.json
This is a read-only check. It regenerates all 360 values at 600 bits, checks
containment in the stored decimal intervals, and independently tests the
integral at three arguments using elementary ball quadrature at 480 bits.
It does not replace the table's original generator or its full certification.
"""
import json
import sys
from decimal import Decimal
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION
def decimal_interval(text):
value = Decimal(text)
midpoint = QQ(Fraction(value))
ulp = QQ(10) ** int(value.as_tuple().exponent)
assert len(value.as_tuple().digits) == 100, (text, value.as_tuple())
return midpoint, ulp
def regenerated_values(doc, bits=600):
C = ComplexBallField(bits)
i = C(0, 1)
rotation = (1 + i) / C(2).sqrt()
for order, rows in doc['Numbers'].items():
n = int(order)
for argument, text in rows.items():
value = ((-i)**n * (C(QQ(argument))*rotation).bessel_K(n)).real()
yield order, argument, text, value
def main():
doc = json.load(sys.stdin)
assert list(doc['Parameters']) == ['nu', 'x']
assert list(doc['Numbers']) == ['0', '1', '2']
expected = sorted({QQ(a)/b for b in range(1, 5) for a in range(1, 20*b+1)})
assert len(expected) == 120
for rows in doc['Numbers'].values():
assert [QQ(x) for x in rows] == expected
max_radius = 0
max_ulp_error = 0
count = 0
for n, x, text, ball in regenerated_values(doc):
midpoint, ulp = decimal_interval(text)
R = ball.parent()
assert ball.is_finite(), (n, x, 'non-finite')
assert R(ball.rad()) < R(ulp)/1000000, (n, x, 'wide ball')
error = abs(ball - R(midpoint)) / R(ulp)
assert error.upper() < 1, (n, x, 'outside written interval')
max_ulp_error = max(max_ulp_error, float(error.upper()))
max_radius = max(max_radius, float(ball.rad()))
count += 1
assert count == 360
print(json.dumps({'regenerated_rows': count, 'precision_bits': 600,
'max_absolute_radius': max_radius,
'max_error_in_written_ulps': max_ulp_error,
'all_rows_have_100_significant_digits': True}), flush=True)
C = ComplexBallField(480)
R = RealBallField(480)
tolerance = R(2)**(-440)
controls = []
for f, start, end, expected in [
(lambda t, analytic: t, C(0), C(1), QQ(1)/2),
(lambda t, analytic: t.sin(), C(0), C.pi()/2, QQ(1)),
]:
answer = C.integral(f, start, end, abs_tol=tolerance, rel_tol=tolerance)
assert answer.real().is_finite() and answer.imag().is_finite() and (answer.real()-R(expected)).contains_zero() and answer.imag().contains_zero()
assert R(answer.real().rad()) < R(2)**(-400)
controls.append(str(answer))
b = R(1)/(8*R(2).sqrt())
u0 = R(9).exp()
tail = (-b*u0).exp()*(u0/b + 1/b**2)
assert tail.upper() < R(QQ(819)/100 * QQ(10)**(-307)).lower()
integral_checks = []
for n, x in [(0, QQ(1)/4), (1, QQ(1)), (2, QQ(20))]:
a = C(x)/C(2).sqrt()
phase = n*C.pi()/2
def integrand(t, analytic):
z = a*t.cosh()
return (-z).exp()*(z+phase).cos()*(n*t).cosh()
value = C.integral(integrand, 0, 9, abs_tol=tolerance, rel_tol=tolerance)
assert value.real().is_finite() and value.imag().is_finite() and value.imag().contains_zero()
midpoint, ulp = decimal_interval(doc['Numbers'][str(n)][str(x)])
error = abs(value.real()-R(midpoint)) + tail
assert R(value.real().rad()) < R(ulp)/1000000
assert error.upper() < R(ulp).lower(), (n,x,'integral mismatch')
integral_checks.append({'nu':n, 'x':str(x), 'radius':str(value.real().rad()),
'contained_with_tail':True})
i = C(0,1)
half = QQ(1)/2
value = ((-i*C.pi()*C(half)/2).exp()*((1+i)/C(2).sqrt()).bessel_K(half)).real()
elementary = (R.pi()/2).sqrt()*(-1/R(2).sqrt()).exp()*(1/R(2).sqrt()+3*R.pi()/8).cos()
assert value.is_finite() and elementary.is_finite()
assert R(value.rad()) < R(2)**(-400) and R(elementary.rad()) < R(2)**(-400)
assert value.overlaps(elementary)
print(json.dumps({'integral_checks':integral_checks, 'control_integrals':controls,
'tail_upper':str(tail.upper()), 'half_order_control':True,
'sage_version':SAGE_VERSION,
'numberdb_version':getattr(numberdb,'__version__','unknown')}),flush=True)
if __name__ == '__main__':
main()