back to table · edit · history · where entries came from · files · download
8292 bytes, as of the version from 2026-10-26 08:47. Recorded here, not run.
"""Independent checks for Kelvin kei -- numberdb.org/T509.
Download generate.py and check.py alongside this file, then run:
sage -python checks.py
sage -python checks.py --stored # NUMBERDB_API_KEY required for a draft
Every row is checked against certified quadrature of an elementary-function
integrand (DLMF 10.32.9), not a Bessel or Kelvin function implementation.
The tail proof, controls, all decimal intervals, and all independent ball
centres and radii are printed as JSON lines for reproducible evidence.
Hankel, recurrence and negative-order reflection checks supplement this.
No command in this file writes to NumberDB.
"""
import json
import os
import sys
import time
from decimal import Decimal
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
from check import exactness, measure, names_its_rings
from generate import KelvinKei
C = ComplexBallField(480)
R = RealBallField(480)
D = ComplexBallField(600)
S = RealBallField(600)
TOL = R(2)**-440
b = R(1)/(8*R(2).sqrt())
u0 = R(9).exp()
TAIL = (-b*u0).exp()*(u0/b+1/b**2)
# For real t >= 9, x >= 1/4 and 0 <= nu <= 2, the absolute
# integrand is <= exp(-b exp(t) + 2t). Substitution u=exp(t)
# gives integral_{exp(9)}^infinity u exp(-bu) du, exactly TAIL.
def emit(kind, **data):
print(json.dumps({"kind": kind, **data}), flush=True)
def finite_narrow(v, relative=None, absolute=None):
assert v.is_finite(), str(v)
if relative is not None:
assert QQ(v.rad()) < relative*abs(QQ(v.center())), str(v)
if absolute is not None:
assert QQ(v.rad()) < absolute, str(v)
def decimal_interval(text):
dec = Decimal(text)
tup = dec.as_tuple()
integer = int("".join(map(str,tup.digits))) * (-1 if tup.sign else 1)
ulp = QQ(10)**int(tup.exponent)
center = QQ(integer)*ulp
return center-ulp, center+ulp
def encloses(text, v):
lower, upper = decimal_interval(text)
return lower <= QQ(v.center())-QQ(v.rad()) and QQ(v.center())+QQ(v.rad()) <= upper
def written_ball(text):
lower, upper = decimal_interval(text)
return S((lower+upper)/2).add_error(S((upper-lower)/2).upper())
def pair(nu, x):
z=D(x)*D(1,1)/D(2).sqrt()
return (-(D(0,1)*D(nu)*D.pi()/2)).exp()*z.bessel_K(D(nu))
def integral(nu, x):
assert 0 <= nu <= 2 and x >= QQ(1)/4
a=C(x)/C(2).sqrt()
phase=C(nu)*C.pi()/2
def integrand(t, analytic):
# exp, sin and cosh are entire, so this is analytic on every ball.
h=t.cosh()
return -(-a*h).exp()*(a*h+phase).sin()*(C(nu)*t).cosh()
v=C.integral(integrand,0,9,abs_tol=TOL,rel_tol=TOL)
assert v.imag().is_finite() and v.imag().contains_zero()
result=v.real().add_error(TAIL.upper())
finite_narrow(result, absolute=QQ(10)**-125)
return result
def controls():
for label,f,end,expected in (
("integral t from 0 to 1",lambda t,_:t,1,R(1)/2),
("integral sin from 0 to pi/2",lambda t,_:t.sin(),C.pi()/2,R(1))):
v=C.integral(f,0,end,abs_tol=TOL,rel_tol=TOL)
finite_narrow(v.real(),absolute=QQ(10)**-125)
assert v.real().overlaps(expected) and v.imag().contains_zero()
emit("control",name=label,value=str(v),radius=str(v.real().rad()))
half=integral(QQ(1)/2,QQ(1))
# DLMF 10.39.2 after multiplying by exp(-pi*i/4):
# kei_(1/2)(1)=-sqrt(pi/2)*exp(-1/sqrt(2))*sin(1/sqrt(2)+3pi/8).
known=-(R.pi()/2).sqrt()*(-1/R(2).sqrt()).exp()*(1/R(2).sqrt()+3*R.pi()/8).sin()
finite_narrow(known,absolute=QQ(10)**-125)
assert half.overlaps(known)
emit("control",name="elementary half-order kei at x=1",value=str(half),expected=str(known))
# Negative control makes the agreement test demonstrably discriminating.
assert not half.overlaps(known+QQ(1)/10**100)
assert not encloses("0.500000000000000000000000000000",R(3)/4)
emit("negative-controls",passed=2)
def records(document):
out={}
for n,group in document["Numbers"].items():
for x,entry in group.items():
out[(str(n),str(x))]=entry["number"] if isinstance(entry,dict) else entry
return out
def run(stored=False):
start=time.monotonic()
emit("versions",sage=sage_version,numberdb=numberdb.__version__,python=sys.version.split()[0])
emit("tail-bound",bound=str(TAIL),cutoff=9,b=str(b))
controls()
generator=KelvinKei()
params=list(generator.enumerate())
values={(p["nu"],p["x"]):generator.value(p,100) for p in params}
assert len(values)==360
assert not names_its_rings(os.path.join(os.path.dirname(__file__),"generate.py"))
assert not exactness(values)
generated={k:to_text(v,100) for k,v in values.items()}
texts=records(numberdb.table("T509")) if stored else generated
assert set(texts)==set(generated)
# Use the actual sibling values in the recurrence, not our own ker code.
ker=records(numberdb.table("T493"))
max_radius=max(v.rad() for v in values.values())
max_relative=max(S(v.rad())/abs(S(v.center())) for v in values.values())
max_integral_radius=R(0)
max_hankel_radius=R(0)
for index,(key,text) in enumerate(texts.items(),1):
ns,xs=key; n=int(ns); x=QQ(xs); value=values[key]
finite_narrow(value,relative=QQ(10)**-100)
assert encloses(text,value),("generator does not fit written interval",key,text,str(value))
independently=integral(QQ(n),x)
assert encloses(text,independently),("integral does not fit written interval",key,text,str(independently))
assert independently.overlaps(value),("generator and integral disjoint",key)
max_integral_radius=max(max_integral_radius,R(independently.rad()))
w=D(x)*D(-1,1)/D(2).sqrt()
hankel=(D.pi()*D(0,1)/2*(w.bessel_J(D(n))+D(0,1)*w.bessel_Y(D(n)))).imag()
finite_narrow(hankel,absolute=QQ(10)**-125)
assert encloses(text,hankel),("Hankel",key)
max_hankel_radius=max(max_hankel_radius,R(hankel.rad()))
negative=pair(-n,x).imag()*(-1)**n
finite_narrow(negative,absolute=QQ(10)**-125)
assert encloses(text,negative),("reflection",key)
center=written_ball(text)
def neighboring(m):
if (str(m),xs) in texts:
return written_ball(texts[(str(m),xs)])
v=pair(m,x).imag()
finite_narrow(v,absolute=QQ(10)**-125)
return v
left=neighboring(n-1)+neighboring(n+1)
right=-S(n)*S(2).sqrt()/S(x)*(center+written_ball(ker[key]))
finite_narrow(left,absolute=QQ(10)**-90)
finite_narrow(right,absolute=QQ(10)**-90)
assert (left-right).contains_zero(),("recurrence",key)
emit("entry",nu=ns,x=xs,number=text,
generator_center=str(value.center()),generator_radius=str(value.rad()),
integral_center=str(independently.center()),integral_radius=str(independently.rad()),
hankel_radius=str(hankel.rad()),passed=["containment","quadrature","hankel","reflection","recurrence"])
# The order-zero origin limit follows analytically from DLMF 10.65.4:
# bei_0(x)=O(x^2), ber_0(x)=1+O(x^4), and the remaining series is O(x^2).
# The following check supplements (rather than proves) that limit.
for k in (20,40,80):
x=QQ(2)**-k
residual=pair(0,x).imag()+S.pi()/4
finite_narrow(residual,absolute=QQ(10)**-125)
assert abs(residual) < S(2)**(-2*k)*k
emit("summary",mode="stored" if stored else "computed",entries=len(texts),
controls=3,negative_controls=2,integral_checks=360,hankel_checks=360,
reflection_checks=360,recurrence_checks=360,origin_limit_checks=3,
max_generator_radius=str(max_radius),max_relative_radius=str(max_relative),
max_integral_radius=str(max_integral_radius),max_hankel_radius=str(max_hankel_radius),
written_size=measure(texts),seconds=round(time.monotonic()-start,2),passed=True)
if __name__=="__main__":
from_stdin=os.environ.get("NUMBERDB_KEY_FROM_STDIN")=="1"
if from_stdin:
os.environ["NUMBERDB_API_KEY"]=sys.stdin.read().strip()
run(stored=from_stdin or "--stored" in sys.argv)