back to table · edit · history · where entries came from · files · download
3118 bytes, as of the version from 2026-10-27 01:14 (current). Recorded here, not run.
"""Check T509 decimals independently of its K-based generator.
Run: sage -python check_t509_repair.py T509-before.json
When no path is supplied, T509-before.json is read beside this file.
No network access or API key is required. The 30 comparisons use the
Hankel J/Y representation at 600 bits and require finite narrow balls
contained in the stored decimal intervals. This is a spot check, not a
rerun of the builder's 360-entry quadrature certificate.
"""
import json
import sys
from pathlib import Path
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 rational(s):
f = Fraction(str(s))
return QQ(f.numerator) / QQ(f.denominator)
def interval_check(value, written):
centre = rational(written)
ulp = QQ(10) ** int(Decimal(written).as_tuple().exponent)
assert value.is_finite()
assert R(value.rad()) < R(ulp)/1000
assert abs(value - R(centre)) < R(ulp)
assert not abs(value - R(centre + 3*ulp)) < R(ulp)
return float(R(abs(value - R(centre)).upper())/R(ulp))
path = Path(sys.argv[1]) if len(sys.argv)>1 else Path(__file__).with_name('T509-before.json')
doc = json.loads(path.read_text())
R = RealBallField(600)
C = ComplexBallField(600)
i = C(0, 1)
pi = C.pi()
rows = doc['Numbers']
assert sum(len(v) for v in rows.values()) == 360
expected = {str(QQ(a)/QQ(b)) for b in range(1,5) for a in range(1,20*b+1)}
assert len(expected) == 120
assert set(rows) == {'0','1','2'}
assert all(set(r)==expected for r in rows.values())
# An elementary half-order control tests phase and rotation before the rows.
x = C(1)
z = x*(-1+i)/C(2).sqrt()
h = (pi*i/2 * (z.bessel_J(QQ(1)/2) + i*z.bessel_Y(QQ(1)/2))).imag()
elementary = -(R.pi()/2).sqrt() * (-1/R(2).sqrt()).exp() * (1/R(2).sqrt()+3*R.pi()/8).sin()
assert h.is_finite() and elementary.is_finite()
assert R(h.rad()) < R('1e-170') and R(elementary.rad()) < R('1e-170')
assert (h-elementary).contains_zero()
checks=[]
for nu in range(3):
for xstr in ['1/4','1/3','1/2','1','2','3','4','5','10','20']:
z=C(rational(xstr))*(-1+i)/C(2).sqrt()
v=(pi*i/2*(z.bessel_J(nu)+i*z.bessel_Y(nu))).imag()
ratio=interval_check(v,rows[str(nu)][xstr])
checks.append({'nu':nu,'x':xstr,'error_upper_in_ulps':ratio})
# Execute the actual displayed Sage program, then compare its returned ball.
namespace={}
exec(doc['Programs']['program-sage']['code'].replace('print(((-i*C.pi()*nu/2).exp()*z.bessel_K(nu)).imag())','result = ((-i*C.pi()*nu/2).exp()*z.bessel_K(nu)).imag()'),namespace)
program_ratio=interval_check(R(namespace['result']),rows['0']['1'])
print(json.dumps({'sage_version':SAGE_VERSION,'precision_bits':600,'entries':360,'arguments_per_order':120,'hankel_checks':checks,'maximum_error_upper_in_ulps':max(x['error_upper_in_ulps'] for x in checks),'half_order_control':True,'shifted_value_controls_rejected':30,'displayed_sage_program_error_upper_in_ulps':program_ratio},indent=2))