back to table · edit · history · where entries came from · files · download
7107 bytes, as of the version from 2026-10-28 19:40 (current). Recorded here, not run.
"""Independent numerical checks of T588, positive ker_0 zeros.
sage -python checks.py # read draft/table and check stored values
sage -python checks.py --local # check freshly generated written values
Set NUMBERDB_API_KEY for a draft. The independent real-series calculation
uses no Bessel or Kelvin calls and explicit geometric tails (root-counting.md).
The additional Arb Bessel calculations check identities and controls only.
"""
import json
import os
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.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_mpfi import RealIntervalField
import sage.version
import mpmath
SERIES_BITS = 1024
R = RealBallField(SERIES_BITS)
I = RealIntervalField(SERIES_BITS)
C = ComplexBallField(SERIES_BITS)
S = C(1,1)/R(2).sqrt()
def qdecimal(text):
f = Fraction(str(text))
return QQ(f.numerator)/f.denominator
def strict_sign(x):
assert x.is_finite(), 'nonfinite result'
if x > 0:
return 1
if x < 0:
return -1
raise AssertionError('zero lies in an endpoint enclosure')
def independent_series(x, threshold=QQ(1)/10**220):
"""DLMF 10.65.2/5 as real arithmetic with an explicit tail bound."""
assert x > 0
t = R(x*x/4)
L = -(R(x)/2).log()-R.euler_constant()
term = R(1)
h = R(0)
ber, bei, hr, hi = R(1), R(0), R(0), R(0)
for k in range(1,10001):
term = term*t/(k*k)
h += R(1)/k
if k%4 == 0:
ber += term
hr += term*h
elif k%4 == 1:
bei += term
hi += term*h
elif k%4 == 2:
ber -= term
hr -= term*h
else:
bei -= term
hi -= term*h
next_term = term*t/((k+1)*(k+1))
ratio = t/((k+2)*(k+2))
if ratio < QQ(1)/2:
plain_tail = next_term/(1-ratio)
harmonic_tail = next_term*((h+R(1)/(k+1))/(1-ratio)
+ratio/((k+2)*(1-ratio)**2))
error = (abs(L)+R.pi()/4)*plain_tail+harmonic_tail
if error < threshold:
break
else:
raise AssertionError('series did not reach its tail target')
error_ball = R(I(-error.upper(),error.upper()))
real = L*ber+R.pi()*bei/4+hr+error_ball
imag = L*bei-R.pi()*ber/4+hi+error_ball
for v in (real,imag):
assert v.is_finite() and v.rad() < QQ(1)/10**200
return real, imag, k
def decimal_interval(text):
d = Decimal(text)
exponent = d.as_tuple().exponent
ulp = QQ(10)**int(exponent)
center = qdecimal(text)
assert len(d.as_tuple().digits) == 100
return center, ulp
def controls():
# A nonzero exact closed form tests the evaluator before any agreement.
known = (R.pi()/2).sqrt()*R(-1).exp()
measured = C(1).bessel_K(QQ(1)/2).real()
assert known > QQ(4)/10 and known < QQ(5)/10
assert measured > 0 and measured.overlaps(known)
assert measured.rad() < QQ(1)/10**280
# A sign known analytically from the origin exclusion bound.
sr,si,_ = independent_series(QQ(1)/4)
assert sr > 1 and si < 0
direct = (S/4).bessel_K(0)
assert sr.overlaps(direct.real()) and si.overlaps(direct.imag())
# Origin expansion with a quantitative error bound, not just tiny digits.
x = QQ(1)/2**20
sr,si,_ = independent_series(x)
L = -(R(x)/2).log()-R.euler_constant()
t = x*x/4
bound = (L+R.pi()/4)*t/(1-t)+t/(1-t)**2
assert abs(sr-L) < bound and abs(si+R.pi()/4) < bound
return {'closed_form_K_half_at_one':str(measured),
'origin_positive_control':True,'origin_error_bound_control':True}
def check_records(records):
evidence = {'table':'T588','instruction_manifest':
'agents/runs/20261028T190129Z-build-skills/manifest.json',
'skill_sha256_local':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
'skill_sha256_live':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
'versions':{'Sage':sage.version.version,'numberdb':numberdb.__version__,
'mpmath':mpmath.__version__},
'series_bits':SERIES_BITS,'controls':controls(),'entries':[]}
assert [str(row['n']) for row in records] == [str(n) for n in range(1,21)]
previous = QQ(0)
controls_as = ['1.71854','6.12728','10.56294','15.00269','19.44381']
for row in records:
n = int(row['n']); text = row['number']
center, ulp = decimal_interval(text)
assert center-ulp > previous
previous = center+ulp
left, _, nl = independent_series(center-ulp)
right, _, nr = independent_series(center+ulp)
sl, sr = strict_sign(left), strict_sign(right)
assert sl*sr == -1
assert sl == (-1)**(n-1)
# Ten last-place units away must NOT pass as a root interval.
displaced_left,_,_ = independent_series(center+9*ulp)
displaced_right,_,_ = independent_series(center+11*ulp)
assert strict_sign(displaced_left)==strict_sign(displaced_right)
# Analytic identities used in the count, checked at every listed root.
z = S*center
g = z.bessel_K(0)
gp = -S*z.bessel_K(1)
gpp = S*S*(z.bessel_K(0)+z.bessel_K(2))/2
residual = center**2*gpp+center*gp-C(0,1)*center**2*g
assert residual.real().contains_zero() and residual.imag().contains_zero()
assert residual.real().rad() < QQ(1)/10**180
assert residual.imag().rad() < QQ(1)/10**180
w = g.real()*gp.imag()-g.imag()*gp.real()
assert w < 0
if n <= 5:
assert abs(center-qdecimal(controls_as[n-1])) < QQ(1)/10**5
evidence['entries'].append({'n':n,'number':text,'left_sign':sl,'right_sign':sr,
'left_series_radius':str(left.rad()),'right_series_radius':str(right.rad()),
'series_terms':max(nl,nr),'displaced_interval_rejected':True,
'ode_checked':True,'W_negative':True,'decimal_ulp':str(ulp)})
print('Checked n=%d: written bracket, negative control, ODE and phase sign' % n,flush=True)
evidence['passed'] = True
print('EVIDENCE_JSON '+json.dumps(evidence),flush=True)
return evidence
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
raw=sys.stdin.read().strip()
if '=' in raw and raw.split('=',1)[0].strip().isidentifier():
raw=raw.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=raw)
if '--local' in sys.argv:
from generate import KelvinKerZeros
from numberdb._write import to_text
gen=KelvinKerZeros()
rows=[{'n':p['n'],'number':to_text(gen.value(p,100),100)} for p in gen.enumerate()]
else:
doc=numberdb.table('T588')
rows=[{'n':n,'number':(row['number'] if isinstance(row,dict) else row)} for n,row in doc['Numbers'].items()]
check_records(rows)