back to table · edit · history · where entries came from · files · download
7591 bytes, as of the version from 2026-10-27 06:00 (current). Recorded here, not run.
"""Independent numerical checks for numberdb.org/T539.
Run: sage -python checks.py, beside generate.py, check.py and dry_run.py.
Those last two files are the campaign's unmodified common checks, attached for
reproducibility. Offline runs need no key. Set NUMBERDB_KEY_FROM_STDIN=1 and
pipe an authorized key to also check every stored interval and call verify().
The independent evaluation uses the ODE y''=x^N y, not the determinant's
Gamma formula. With m=N+2, the initial-value solution b(0)=0,b'(0)=1 is
b(x)=sum(c_k*x^(mk+1)), c_0=1, c_k=c_(k-1)/(mk(mk+1)).
At x=1 both series are positive. The ratios for b and b' decrease and are
at most 1/2 for m>=2. After term k, twice the first omitted term bounds each
remaining tail. Both tails are enclosed with exact rational endpoints.
Voros (2000), eqs. (47),(135), gives the normalized decaying solution
psi(x)=2*sqrt(nu/pi)*sqrt(x)*K_nu(2*nu*x^(1/(2*nu))).
Since W(psi,b) is constant, psi(0)=psi(1)*b'(1)-psi'(1)*b(1).
K_nu'(z)=-(K_(nu-1)(z)+K_(nu+1)(z))/2. This computes the determinant from
Bessel values away from the origin plus a rational Taylor solution. It shares
the published spectral normalization with the generator, but no Gamma formula
or implementation code. It does not claim to reprove the spectral theorem.
The m=2 control uses psi=exp(-x), b=sinh(x), so the boundary value must be 1.
The N=1 control uses 2*sqrt(pi)*Ai(0); N=2 additionally differentiates the
independent harmonic spectral zeta 4^(-s)*zeta(s,3/4), using the explicit
Hurwitz derivative argument at 150 and 200 decimal digits in mpmath. The latter is agreement evidence, not an enclosure.
"""
import json
import os
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.real_mpfi import RealIntervalField
from sage.rings.complex_arb import ComplexBallField
import mpmath as mp
from generate import DirichletDeterminants, authenticate_from_stdin
from check import exactness, measure, agrees_with, names_its_rings
import dry_run
def narrow(ball, bits=360):
assert ball.is_finite(), 'Nonfinite ball'
assert ball.accuracy() >= bits, 'Insufficient relative accuracy'
return ball
def independent_boundary(degree, bits=620):
field = RealBallField(bits)
interval = RealIntervalField(bits)
m = degree + 2
nu = field(QQ(1)/m)
c = QQ(1)
b, bp = QQ(1), QQ(1)
k = 0
threshold = QQ(1) / 10**165
while True:
j = k + 1
omitted = c / QQ((m*j)*(m*j+1))
omitted_derivative = (m*j+1)*omitted
if 2*omitted_derivative < threshold:
break
c = omitted
b += c
bp += (m*j+1)*c
k = j
b = field(interval(b, b+2*omitted))
bp = field(interval(bp, bp+2*omitted_derivative))
cfield = ComplexBallField(bits)
z = cfield(2*nu)
kval_complex = z.bessel_K(cfield(nu))
kd_complex = -(z.bessel_K(cfield(nu-1))+z.bessel_K(cfield(nu+1)))/2
assert kval_complex.imag().contains_zero() and kd_complex.imag().contains_zero()
kval, kd = kval_complex.real(), kd_complex.real()
prefactor = 2*(nu/field.pi()).sqrt()
psi = prefactor*kval
psip = prefactor*(kval/2+kd)
return narrow(psi*bp-psip*b, 520), k+1
def decimal_interval(text, bits=700):
# Exact decimal parsing; Sage QQ('1.25') is not portable with named imports.
from decimal import Decimal
from fractions import Fraction
d = Decimal(text)
f = Fraction(d)
centre = QQ(f.numerator)/f.denominator
exponent = d.as_tuple().exponent
unit = QQ(10)**exponent
return RealIntervalField(bits)(centre-unit, centre+unit)
def run():
authenticate_from_stdin()
path = Path(__file__).with_name('generate.py')
assert dry_run.main([str(path)]) == 0
assert not names_its_rings(str(path))
gen = DirichletDeterminants()
values = {int(p['N']): gen.value(p,100) for p in gen.enumerate()}
assert list(values) == list(range(1,13))
assert not exactness(values)
field = RealBallField(620)
control, _ = independent_boundary(0)
assert 1 in control
assert 0 not in control
assert 2 not in control
print('CONTROL: m=2 boundary is 1, excluding 0 and 2.')
other = {degree: independent_boundary(degree) for degree in values}
assert not agrees_with(values, lambda degree: other[degree][0])
for degree,value in values.items():
narrow(value,400)
assert value > 0
assert value.overlaps(narrow(gen.value({'N':degree},200),700))
assert value.overlaps(other[degree][0])
assert not (value+QQ(1)/10**90).overlaps(other[degree][0])
airy = narrow(2*field.pi().sqrt()*ComplexBallField(620)(0).airy_ai().real(),600)
assert values[1].overlaps(airy)
harmonic = narrow(field.pi().sqrt()/field(QQ(3)/4).gamma(),600)
assert values[2].overlaps(harmonic)
harmonic_zeta = []
for dps in [150,200]:
with mp.workdps(dps):
a = mp.mpf(3)/4
# derivative=1 avoids finite-difference cancellation across the
# Hurwitz reflection branch at negative infinitesimal s.
value = mp.exp(mp.log(4)*mp.zeta(0,a)-mp.zeta(0,a,derivative=1))
harmonic_zeta.append(mp.nstr(value,110))
discrepancy = abs(value-mp.mpf(str(harmonic.mid())))
print('HARMONIC_ZETA_CHECK',dps,'decimal digits; discrepancy',mp.nstr(discrepancy,8))
assert discrepancy < mp.mpf('1e-145')
assert harmonic_zeta[0] == harmonic_zeta[1]
from numberdb._write import to_text
encoded = {degree:to_text(value,digits=100) for degree,value in values.items()}
for degree,text in encoded.items():
interval = decimal_interval(text)
assert RealIntervalField(700)(other[degree][0]) in interval
summary = {'table':'T539', 'entries':12,
'working_bits':numberdb.bits(100,losing=96),
'minimum_generator_accuracy_bits':min(int(v.accuracy()) for v in values.values()),
'minimum_independent_accuracy_bits':min(int(v[0].accuracy()) for v in other.values()),
'series_term_counts':{str(k):v[1] for k,v in other.items()},
'raw_ball_measurement':measure(values),
'stored_decimal_measurement':measure(encoded),
'harmonic_zeta_110_digits':harmonic_zeta[0],
'all_controls_and_perturbation_checks_passed':True,
'all_rounding_intervals_contain_independent_enclosures':True,
'software':{**gen.environment(),'mpmath':mp.__version__},
'values':encoded}
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
report = gen.verify(sample=None)
print(report)
assert report.ok
from numberdb._generate import _stored_entries
observed = _stored_entries(numberdb.table('T539'))
assert len(observed) == 12
# The shared stored() requires Django; use the supported public API here.
for identity, text in observed.items():
degree = int(identity)
interval = decimal_interval(str(text))
assert RealIntervalField(700)(other[degree][0]) in interval
assert RealIntervalField(700)(gen.value({'N':degree},200)) in interval
if degree == 1: assert RealIntervalField(700)(airy) in interval
if degree == 2: assert RealIntervalField(700)(harmonic) in interval
summary['stored_checks_passed'] = True
print('EVIDENCE_JSON='+json.dumps(summary,sort_keys=True))
if __name__ == '__main__':
run()