back to table · edit · history · where entries came from · files · download
8060 bytes, as of the version from 2026-10-25 16:01. Recorded here, not run.
"""Independent checks for Kelvin ker values -- numberdb.org/T493.
sage -python checks.py
sage -python checks.py --stored # NUMBERDB_API_KEY required for the draft
The certified check integrates elementary functions, without calling a
Bessel or Kelvin implementation. See DLMF 10.32.9. At a=x/sqrt(2),
ker_n(x) is the integral over t>=0 of
exp(-a*cosh(t))*cos(a*cosh(t)+n*pi/2)*cosh(n*t).
For 0<=n<=2 and t>=9, the absolute integrand is bounded by
exp(-b*exp(t)+2*t), where b=1/(8*sqrt(2)) since x>=1/4.
Substituting u=exp(t) gives the uniform omitted-tail bound
exp(-b*exp(9))*(exp(9)/b+1/b**2) < 8.19e-307.
The integral on [0,9] is enclosed by 480-bit Arb quadrature. The integrand
is entire, so the analytic flag introduces no branch or endpoint issue.
All comparisons require finite, narrow balls; an overlap alone is never
accepted as evidence. Written decimals are interpreted as plus or minus one
unit in their last place, with exact rational endpoints. Every such interval
must strictly contain both the generator enclosure and the integral enclosure.
mpmath at 160 digits is a supplementary independent-library check, not the
justification for certification. Hankel rotation, reflection and recurrence
checks run on every applicable row; stored mode uses database decimals.
"""
import json
import os
import sys
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from sage.version import version as sage_version
from numberdb._write import to_text
import mpmath as mp
from generate import KelvinKer, arguments, kelvin_pair
from check import exactness, measure, names_its_rings
C = ComplexBallField(480)
R = RealBallField(600)
H = ComplexBallField(600)
B = C(0).real().parent()
TOL = B(2)**-440
b = B(1)/(8*B(2).sqrt())
u0 = B(9).exp()
TAIL = (-b*u0).exp()*(u0/b + 1/b**2)
def narrow(ball, relative=False):
assert ball.is_finite(), ('nonfinite', ball)
assert R(ball.rad()) < R(10)**-125, ('wide', ball)
if relative:
assert R(ball.rad()) < abs(R(ball))*R(10)**-110, ('relative width', ball)
def integral(nu, x):
a = C(QQ(x))/C(2).sqrt()
order = C(QQ(nu))
angle = order*C.pi()/2
def integrand(t, analytic):
ch = t.cosh()
return (-a*ch).exp()*(a*ch+angle).cos()*(order*t).cosh()
result = C.integral(integrand, 0, 9, abs_tol=TOL, rel_tol=TOL)
assert result.imag().is_finite() and result.imag().contains_zero()
value = result.real().add_error(TAIL)
narrow(value, relative=True)
return value
def endpoints(text):
dec = Decimal(text)
center = QQ(Fraction(dec))
unit = QQ(Fraction(10)**dec.as_tuple().exponent)
assert len(dec.as_tuple().digits) == 100, ('not 100 significant digits', text)
return center-unit, center+unit, center, unit
def contains(text, value):
lo, hi, _, _ = endpoints(text)
assert value.is_finite()
assert R(lo) < R(value) and R(value) < R(hi), (text, value)
def as_ball(text):
_, _, center, unit = endpoints(text)
return R(center).add_error(R(unit))
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
stored_mode = '--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED') == '1'
generator = KelvinKer()
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
assert TAIL < B('8.19e-307')
# Controls precede every use of the quadrature as a measurement.
linear = C.integral(lambda t, analytic:t, 0, 1).real()
sine = C.integral(lambda t, analytic:t.sin(), 0, C.pi()/2).real()
assert (QQ(1)/2 in linear) and (1 in sine)
narrow(linear)
narrow(sine)
a = R(1)/R(2).sqrt()
half_expected = (R.pi()/2).sqrt()*(-a).exp()*(a+3*R.pi()/8).cos()
half = integral(QQ(1)/2, 1)
assert half.overlaps(half_expected) and not half.contains_zero()
print('CONTROLS PASS: integral(t)=1/2; integral(sin)=1; elementary half-order ker(1)',flush=True)
mp.mp.dps=160
computed = {}
written = {}
if stored_mode:
doc = numberdb.table(generator.table)
assert list(doc['Parameters']) == ['nu','x']
for n, entries in doc['Numbers'].items():
for x, entry in entries.items():
written[(str(n),str(x))] = entry['number'] if isinstance(entry,dict) else entry
params_list = list(generator.enumerate())
expected_keys={(p['nu'],p['x']) for p in params_list}
if stored_mode:
assert set(written)==expected_keys
evidence=[]
max_radius=R(0)
max_relative_radius=R(0)
max_integral_radius=R(0)
max_mpmath_error=mp.mpf(0)
for j, params in enumerate(params_list):
key=(params['nu'],params['x'])
n=int(params['nu']); x=QQ(params['x'])
value=generator.value(params,100)
computed[key]=value
assert value.is_finite() and not value.contains_zero()
radius=R(value.rad())
relative_radius=radius/abs(R(value))
assert relative_radius < R(10)**-100, (key, relative_radius)
max_radius=max(max_radius,radius)
max_relative_radius=max(max_relative_radius,relative_radius)
text=written[key] if stored_mode else to_text(value,100)
written[key]=text
contains(text,value)
reference=integral(n,x)
contains(text,reference)
max_integral_radius=max(max_integral_radius,R(reference.rad()))
mp_x=mp.mpf(int(x.numerator()))/int(x.denominator())
independent=mp.ker(n,mp_x)
lo,hi,center,unit=endpoints(text)
mp_center=mp.mpf(int(center.numerator()))/int(center.denominator())
mp_unit=mp.mpf(int(unit.numerator()))/int(unit.denominator())
error=abs(independent-mp_center)/mp_unit
assert error < 1, ('mpmath disagreement',key,error)
max_mpmath_error=max(max_mpmath_error,error)
# Independent evaluation via J and Y at the other rotated argument.
w=H(x)*(-1+H(0,1))/H(2).sqrt()
hankel=(H.pi()*H(0,1)/2*(w.bessel_J(n)+H(0,1)*w.bessel_Y(n))).real()
narrow(hankel,relative=True)
contains(text,hankel)
negative=kelvin_pair(-n,x,600).real()*((-1)**n)
narrow(negative,relative=True)
contains(text,negative)
evidence.append({'nu':str(n),'x':str(x),'written':text,'generator_ball':str(value),'integral_ball':str(reference),'hankel_ball':str(hankel),'mpmath':mp.nstr(independent,160)})
if (j+1)%60==0:
print('CHECKED %d/360'% (j+1),flush=True)
assert not exactness(computed)
# Use stored decimal intervals wherever a neighboring order is present.
for params in params_list:
n=int(params['nu']); x=QQ(params['x'])
def order_value(m):
key=(str(m),str(x))
return as_ball(written[key]) if key in written else R(kelvin_pair(m,x,600).real())
kei=kelvin_pair(n,x,600).imag()
narrow(kei)
residual=order_value(n-1)+order_value(n+1)+R(n)*R(2).sqrt()/R(x)*(order_value(n)-kei)
assert residual.is_finite() and residual.contains_zero()
assert R(residual.rad()) < R(10)**-90
sizes=measure(written)
summary={'table':'T493','mode':'stored' if stored_mode else 'local','entries':len(written),'orders':3,'arguments':len(arguments()),'max_generator_radius':str(max_radius),'max_relative_radius':str(max_relative_radius),'max_integral_radius':str(max_integral_radius),'uniform_tail_bound':str(TAIL),'max_mpmath_error_in_last_places':mp.nstr(max_mpmath_error,12),'size':sizes,'controls':3,'integral_checks':360,'mpmath_checks':360,'hankel_checks':360,'negative_order_checks':360,'recurrence_checks':360,'versions':{'sage':sage_version,'numberdb':numberdb.__version__,'mpmath':mp.__version__}}
print('SUMMARY_JSON '+json.dumps(summary),flush=True)
print('EVIDENCE_JSON '+json.dumps(evidence),flush=True)
if __name__=='__main__':
main()