back to table · edit · history · where entries came from · files · download
6461 bytes, as of the version from 2026-10-28 20:01. Recorded here, not run.
"""Independent checks of ber_0 positive roots -- numberdb.org/T589.
sage -python checks.py # emit JSON evidence, no API writes
sage -python checks.py table.json # check saved live Numbers instead
Only candidate enclosures and the root-count evidence are imported from the
generator. Endpoint checks use an explicit real series with a geometric tail,
not its Bessel evaluator. All error claims are interval claims.
"""
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.env import SAGE_VERSION
from fractions import Fraction
from decimal import Decimal
from pathlib import Path
import json
import sys
import hashlib
import mpmath as mp
import generate
from numberdb._generate import _written
def q(text):
f=Fraction(text)
return QQ(f.numerator)/QQ(f.denominator)
def real_series(x, bits=900):
R=RealBallField(bits)
x=R(x)
t=R(1);s=t
threshold=R(QQ(1)/10**180)
for m in range(1000):
ratio=x**4/(16*((2*m+1)*(2*m+2))**2)
if ratio < 1:
tail=abs(t)*ratio/(1-ratio)
if tail < threshold:
result=s.add_error(tail)
assert result.is_finite()
return result,m,str(tail)
t=-t*ratio
s+=t
raise ArithmeticError('series tail failed to converge')
def strict_sign(y):
assert y.is_finite()
if y>0:return 1
if y<0:return -1
raise AssertionError('independent series sign is unresolved')
def decimal_interval(text):
d=Decimal(text)
centre=q(text)
ulp=QQ(10)**int(d.as_tuple().exponent)
return centre,ulp
def check_rows(rows):
R=RealBallField(900);C=ComplexBallField(900)
zero,_,_=real_series(QQ(0))
assert zero == 1
one,_,_=real_series(QQ(1))
assert one > QQ(984)/1000 and one < QQ(985)/1000
c=C(1,1)/R(2).sqrt();j=C(-1,1)/R(2).sqrt()
printed=['2.84892','7.23883','11.67396','16.11356','20.55463']
evidence=[];previous=QQ(0)
assert [int(row['params']['n']) for row in rows]==list(range(1,21))
for row in rows:
n=int(row['params']['n']);text=row['number']
centre,ulp=decimal_interval(text)
assert len(Decimal(text).as_tuple().digits)==100
assert centre-ulp > previous
previous=centre+ulp
left,k1,tail1=real_series(centre-ulp)
right,k2,tail2=real_series(centre+ulp)
assert strict_sign(left)==(-1)**(n-1)
assert strict_sign(right)==(-1)**n
for v in [left,right]:
assert R(v.rad()) < abs(v)/10**30
# A shifted false bracket must fail, including at the largest root.
wrong_l,_,_=real_series(centre+9*ulp)
wrong_r,_,_=real_series(centre+11*ulp)
assert strict_sign(wrong_l)==strict_sign(wrong_r)
# Independently evaluated Bessel representations and their derivatives.
fi=(C(centre)*c).bessel_I(0)
fj=(C(centre)*j).bessel_J(0)
dpi=c*(C(centre)*c).bessel_I(1)
dpj=-j*(C(centre)*j).bessel_J(1)
for val in [fi,fj,dpi,dpj]:
assert val.real().is_finite() and val.imag().is_finite()
assert R(val.rad()) < abs(val)*R(QQ(1)/10**200)
assert (fi-fj).contains_zero() and (dpi-dpj).contains_zero()
series_mid,_,_=real_series(centre)
assert (series_mid-fi.real()).contains_zero()
# The ODE and increasing-phase invariant, with J2 derivative recurrence.
fpp=j*j*((C(centre)*j).bessel_J(2)-fj)/2
assert (fpp+dpj/C(centre)-C(0,1)*fj).contains_zero()
W=fj.real()*dpj.imag()-fj.imag()*dpj.real()
assert W>0
if n<=5:
assert abs(centre-q(printed[n-1])) < QQ(5)/10**6
# Independent software and independent precision levels. These are
# agreement controls; the rigorous real-series signs establish digits.
approximations=[]
for dps in [160,210]:
with mp.workdps(dps):
v=mp.findroot(lambda x: mp.hyper([], [mp.mpf(1)/2,mp.mpf(1)/2,1],-x**4/256),
(mp.mpf(text)-mp.mpf('0.01'),mp.mpf(text)+mp.mpf('0.01')))
written=mp.nstr(v,130)
assert abs(q(written)-centre) < ulp
approximations.append(written)
with mp.workdps(130):
seed=mp.sqrt(2)*(n-mp.mpf(3)/8)*mp.pi
program_root=mp.findroot(lambda x: mp.ber(0,x), (seed-mp.mpf('0.25'),seed+mp.mpf('0.25')))
assert abs(q(mp.nstr(program_root,120))-centre)<ulp
evidence.append({'n':n,'number':text,'written_ulp':str(ulp),
'left_sign':strict_sign(left),'right_sign':strict_sign(right),
'left_series':str(left),'right_series':str(right),
'series_terms':max(k1,k2)+1,'tail_bounds':[tail1,tail2],
'displaced_bracket_rejected':True,
'bessel_rotation_and_derivative':True,'ode_and_W_positive':True,
'hypergeometric_roots_160_210':approximations})
return evidence
def main():
if len(sys.argv)>1:
doc=json.loads(Path(sys.argv[1]).read_text())
rows=doc['Numbers']
else:
g=generate.KelvinBerZeros()
rows=[{'params':p,'number':_written(g.value(p,100),100,'decimal')} for p in g.enumerate()]
checks=check_rows(rows)
brackets,cells=generate.counted_brackets()
result={'table':'T589','skill_manifest':'../20261028T195045Z-build-skills/manifest.json',
'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
'versions':{'sage':SAGE_VERSION,'numberdb':numberdb.__version__,'mpmath':mp.__version__},
'working_bits':800,'independent_series_bits':900,'entries':len(rows),
'known_answer_controls':['ber_0(0)=1','0.984<ber_0(1)<0.985',
'five A&S printed roots within 0.000005'],
'counted_cells':len(cells),'brackets':[[str(a),str(b)] for a,b in brackets],
'count_certificate':cells,'entry_checks':checks,
'files':{name:hashlib.sha256(Path(__file__).with_name(name).read_bytes()).hexdigest()
for name in ['generate.py','checks.py','certificate.md','argument-selection.md']},
'all_passed':True}
print(json.dumps(result,indent=2))
if __name__=='__main__':main()