back to table · edit · history · where entries came from · files · download
7275 bytes, as of the version from 2026-10-27 06:27 (current). Recorded here, not run.
"""Independent numerical checks for numberdb.org/T540.
Run: sage -python checks.py
To repeat on the live draft: NUMBERDB_API_KEY=... sage -python checks.py --stored
The campaign wrapper instead uses NUMBERDB_KEY_FROM_STDIN=1 and CHECK_STORED=1.
Requires generate.py and the campaign's unchanged check.py beside this file.
No write endpoint is used.
Independent reconstruction: for m=N+2, let F solve F''=x^N F,
F(0)=1, F'(0)=0. At x=1 its series terms are a_0=1 and
a_k=a_(k-1)/[(km)(km-1)].
The derivative terms are km*a_k. All are positive. After term K the
successive ratios decrease; hence tail <= next_term/(1-next_ratio).
Exact rational sums plus these explicit tails enclose F(1), F'(1).
Voros eqs. (47),(135) give D+ = -psi'(0), with
psi(x)=2*sqrt(nu/pi)*sqrt(x)*K_nu(2*nu*x^((N+2)/2)).
The constant Wronskian gives D+ = F'(1)*psi(1)-F(1)*psi'(1).
This uses Bessel K and an ODE power series, not the generator's Gamma formula.
The N=0 auxiliary control has psi=e^(-x), F=cosh(x), and D+=1.
"""
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.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
import mpmath as mp
from generate import NeumannDeterminants
from check import exactness, measure, names_its_rings
CHECK_BITS = 720
TERMS = 100
def number(v):
return v['number'] if isinstance(v, dict) else v
def relative_bits(v):
assert v.is_finite() and not v.contains_zero()
# Conservative lower bound, using strict ball comparison.
b = 0
while v.parent()(v.rad()) * 2**(b+1) < abs(v):
b += 1
return b
def reconstruction(degree):
B = RealBallField(CHECK_BITS)
C = ComplexBallField(CHECK_BITS)
m = degree + 2
nu = QQ(1)/m
term = QQ(1)
f = QQ(1)
df = QQ(0)
for k in range(1, TERMS+1):
term /= QQ(k*m*(k*m-1))
f += term
df += k*m*term
nxt = term / QQ((TERMS+1)*m*((TERMS+1)*m-1))
ratio = QQ(1) / ((TERMS+2)*m*((TERMS+2)*m-1))
dratio = QQ(1) / ((TERMS+1)*m*((TERMS+2)*m-1))
fball = B(f) + B(0, nxt/(1-ratio))
dfball = B(df) + B(0, (TERMS+1)*m*nxt/(1-dratio))
z = C(2*nu)
vals = [z.bessel_K(order) for order in [nu,nu-1,nu+1]]
assert all(v.real().is_finite() and v.imag().is_finite() and v.imag().contains_zero() for v in vals)
kv, km, kp = [v.real() for v in vals]
prefactor = 2*(B(nu)/B.pi()).sqrt()
psi = prefactor*kv
dpsi = prefactor*(kv-km-kp)/2
return dfball*psi-fball*dpsi
def rational_decimal(s):
q = Fraction(Decimal(s))
return QQ(q.numerator)/q.denominator
def rounding_interval(s):
q = rational_decimal(s)
exponent = Decimal(s).as_tuple().exponent
unit = QQ(10)**exponent
return q-unit, q+unit
def tight(v, minimum=600):
assert v.is_finite() and relative_bits(v) >= minimum
def main():
B = RealBallField(CHECK_BITS)
C = ComplexBallField(CHECK_BITS)
control = reconstruction(0)
tight(control)
assert (1 in control)
assert not control.overlaps(B(1)+B(QQ(1)/10**90))
# Independent Airy boundary and harmonic spectral-zeta anchors.
airy = -2*B.pi().sqrt()*C(0).airy_ai_prime().real()
harmonic = 2*B.pi().sqrt()/B(QQ(1)/4).gamma()
tight(airy); tight(harmonic)
generator = NeumannDeterminants()
values = {p['N']:number(generator.value(p,100)) for p in generator.enumerate()}
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
assert not exactness(values)
report = {'table':'T540','control':{'N':0,'known':1,'relative_bits':relative_bits(control)},
'working_bits':429,'independent_bits':CHECK_BITS,'series_terms':TERMS,
'measure':measure(values),'entries':[]}
independent = {}
for key,v in values.items():
n = int(key)
hi = number(generator.value({'N':key},200))
ind = reconstruction(n)
independent[key]=ind
tight(hi); tight(ind)
assert relative_bits(v)>=400
assert v.overlaps(hi) and v.overlaps(ind)
assert not (v+B(QQ(1)/10**90)).overlaps(ind)
# Product identity with a separately evaluated Dirichlet formula.
nu = QQ(1)/(n+2)
dm = B(nu).gamma()*B(nu)**(QQ(1)/2-nu)/B.pi().sqrt()
identity = dm*ind*(B.pi()*B(nu)).sin()
tight(identity); assert (1 in identity)
if n==1: assert v.overlaps(airy)
if n==2: assert v.overlaps(harmonic)
report['entries'].append({'N':n,'generator_relative_bits':relative_bits(v),
'independent_relative_bits':relative_bits(ind),'value':str(v),
'independent':str(ind),'precision_agrees':True,'ode_agrees':True,
'product_identity':True,'perturbation_rejected':True})
mp_checks=[]
for dps in [150,200]:
with mp.workdps(dps):
# An excessively small default step exposes cancellation in the
# Hurwitz-zeta evaluator near s=0 (see zeta-diagnostic.log).
# Controlled central differences have O(h^2) truncation error;
# this is an extra heuristic check, not the certified computation.
step_digits = 60 if dps == 150 else 80
h = mp.mpf(10)**(-step_digits)
threshold = mp.mpf(10)**(-2*step_digits+5)
derivative = lambda f: mp.diff(f,0,h=h,addprec=100)
assert abs(derivative(lambda s:mp.exp(2*s))-2)<threshold
anchor=mp.exp(-derivative(lambda s:mp.power(4,-s)*mp.zeta(s,mp.mpf(1)/4)))
delta=abs(anchor-mp.mpf(str(harmonic.mid())))
assert delta<threshold
mp_checks.append({'dps':dps,'step':'1e-'+str(step_digits),
'harmonic_zeta_derivative':mp.nstr(anchor,dps),
'difference':mp.nstr(delta,8)})
report['harmonic_zeta_checks']=mp_checks
if '--stored' in sys.argv or os.environ.get('CHECK_STORED')=='1':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
doc=numberdb.table('T540')
assert list(doc['Numbers'])==list(values)
dm_doc=numberdb.table('T539')
for key,entry in doc['Numbers'].items():
s = entry['number'] if isinstance(entry,dict) else entry
lo,hi=rounding_interval(s)
assert independent[key]>B(lo) and independent[key]<B(hi)
# Test the cross-table identity on both actual stored intervals.
dm_s=dm_doc['Numbers'][key]
dm_s=dm_s['number'] if isinstance(dm_s,dict) else dm_s
dl,dh=rounding_interval(dm_s)
lhs=(B((lo+hi)/2, (hi-lo)/2)*B((dl+dh)/2,(dh-dl)/2)
*(B.pi()/ (int(key)+2)).sin())
assert lhs.is_finite() and relative_bits(lhs)>320 and (1 in lhs)
report['stored_rounding_encloses_independent']=True
report['stored_product_identity']=True
verify=generator.verify(sample=None)
assert verify.ok
report['verify']=str(verify)
report['versions']={'numberdb':numberdb.__version__,'mpmath':mp.__version__}
print(json.dumps(report,indent=2))
if __name__=='__main__':
main()