back to table · edit · history · where entries came from · files · download
7880 bytes, as of the version from 2026-10-27 08:39 (current). Recorded here, not run.
"""Independent formula and stored-value checks for numberdb.org/T542.
Run with SageMath: sage -python checks.py [--library] [--stored]
Download generate.py, accelerate.py, check.py and values.json alongside this.
The default compares the generator's saved 200-digit results against a fresh
190-digit evaluation of Voros (1983), (C.22), (C.28), (C.33). The generic
asymptotic summation helper is shared; this is independent formula checking,
not a wholly independent algorithm. --library evaluates the different Voros
series with mpmath's hyper() at 25 digits instead, sharing no tail code.
--stored compares fresh references against the live draft's written intervals
and tests both Wronskian identities reconstructed from stored full and parity
traces. NUMBERDB_API_KEY is needed for private tables.
"""
import os
import sys
import json
from pathlib import Path
import numberdb.sage as numberdb
import mpmath as mp
from sage.rings.real_mpfi import RealIntervalField
from accelerate import hyper_tail
from check import agrees_with, exactness, measure
def voros(N, n, working=190, library=False):
with mp.workdps(working):
if N == 1:
# Log determinant ratio, differentiated independently of the
# Gamma expression used by the generator.
f = lambda E: mp.log(-mp.airyai(-E,1)/mp.airyai(-E))
return mp.nstr(-mp.diff(f, 0, 2),working), None
if N == 2:
# The two arithmetic progressions in the harmonic spectrum.
return mp.nstr((mp.zeta(n,mp.mpf(1)/4)-mp.zeta(n,mp.mpf(3)/4))/4**n,working), None
v = mp.mpf(1)/(N+2); h = mp.mpf(1)/2
P1 = mp.sin(mp.pi*v)/(2*mp.sqrt(mp.pi))*(2*v)**(2-4*v)
P1 *= mp.gamma(v)*mp.gamma(2*v)*mp.gamma(3*v)/mp.gamma(h+2*v)
if n == 1:
return mp.nstr(P1,working), mp.nstr(mp.tan(2*mp.pi*v)/mp.tan(mp.pi*v)*P1,working)
hyper = (lambda A,B: mp.hyper(A,B,1)) if library else (lambda A,B: hyper_tail(A,B,240,160))
F = hyper([4*v,2*v,3*v,5*v,h],[1+2*v,1+v,1-v,h+4*v])
pref = mp.sin(mp.pi*v)/(4*v*v*mp.sqrt(mp.pi))*(2*v)**(4-8*v)
pref *= mp.gamma(3*v)*mp.gamma(4*v)*mp.gamma(5*v)/mp.gamma(h+4*v)
full = P1**2+pref*F
Fp = hyper([6*v,3*v,4*v,5*v,h+v],[1+3*v,1+2*v,1+v,h+5*v])
pref_p = (2*v)**(4-8*v)/(6*v)
pref_p *= mp.gamma(4*v)*mp.gamma(5*v)*mp.gamma(6*v)*mp.gamma(h+v)
pref_p /= mp.gamma(1+v)*mp.gamma(1+2*v)*mp.gamma(h+5*v)
return mp.nstr(full-pref_p*Fp,working),mp.nstr(full,working)
def agreeing_digits(a,b):
a,b=mp.mpf(a),mp.mpf(b)
return -mp.log10(abs((a-b)/b)) if a!=b else mp.inf
def contains_written(text, value):
# Decimal spelling denotes +/- one unit in the last place. This tests a
# high-precision heuristic reference, not an asserted rigorous enclosure.
from decimal import Decimal
t=Decimal(text)
ulp=mp.mpf(10)**int(t.as_tuple().exponent)
return abs(mp.mpf(str(t))-mp.mpf(value))<=ulp
def written_interval(text, field):
from decimal import Decimal
t=Decimal(text)
mid=field(text)
radius=field(10)**int(t.as_tuple().exponent)
return field(mid.lower()-radius.upper(),mid.upper()+radius.upper())
def run_checks(stored_mode=False,library=False):
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if key.startswith('NUMBERDB_') and '=' in key:
key=key.split('=',1)[1].strip().strip(chr(34)+chr(39))
os.environ['NUMBERDB_API_KEY']=key
mp.mp.dps=220
saved=json.loads(Path(__file__).with_name('values.json').read_text())
# Nonzero controls and a negative control precede interpretation.
with mp.workdps(190):
assert abs(hyper_tail([mp.mpf(1),mp.mpf(1)],[mp.mpf(3)],240,160)-2)<mp.mpf('1e-150')
if library:
with mp.workdps(25):
assert abs(mp.hyper([1,1],[3],1)-2)<mp.mpf('1e-24')
assert not contains_written('2.'+'0'*99,mp.mpf(2)+mp.mpf('1e-90'))
# Airy derivative and Hurwitz implementations must return the known
# nonzero harmonic answer, not collapse to zero by a coercion mistake.
assert agreeing_digits(voros(2,2)[0],mp.catalan)>180
print('CONTROLS PASS',flush=True)
from numberdb._generate import _stored_entries
live=_stored_entries(numberdb.table('T542')) if stored_mode else None
full_live=_stored_entries(numberdb.table('T541')) if stored_mode else None
field=RealIntervalField(800)
results={};values={};expected={}
assert len(saved)==22
if live is not None: assert set(live)==set(saved)
for key,row in saved.items():
N,n=map(int,key.split(','))
actual=live[key] if stored_mode else row['value200']
ref,full=voros(N,n,25 if library else 190,library)
d=agreeing_digits(actual,ref)
threshold=20 if library else (98 if stored_mode else 115)
assert d>threshold,(key,'independent formula disagreement',str(d))
assert mp.mpf(actual)>0
if stored_mode and not library:
assert contains_written(actual,ref),(key,'outside stored rounding interval')
residual=None
if N>=3 and not library:
if stored_mode:
# Propagate the actual printed rounding intervals of both
# tables. Recover each parity sector from S and P.
p1=written_interval(live[f'{N},1'],field)
sn=written_interval(full_live[key],field)
pn=written_interval(live[key],field)
zm=(sn-pn)/2;zp=(sn+pn)/2
sigma=field(2)/(N+2)
coeff=lambda a:(field.pi()*sigma*(field(a)+field(1)/2)).sin()/(field.pi()*sigma/2).sin()
if n==1: r=coeff(1)*zm+coeff(-1)*zp
else: r=coeff(2)*zm+coeff(-2)*zp+(coeff(1)**2-coeff(2))*p1**2
assert not r.is_NaN() and r.contains_zero() and r.absolute_diameter()<field('1e-94').lower(),(key,str(r))
residual=str(r)
else:
p1=mp.mpf(saved[f'{N},1']['value200'])
zm=mp.mpf(row['dirichlet200']);zp=mp.mpf(row['neumann200'])
sigma=mp.mpf(2)/(N+2)
coeff=lambda a:mp.sin(mp.pi*sigma*(mp.mpf(a)+mp.mpf(1)/2))/mp.sin(mp.pi*sigma/2)
r=coeff(1)*zm+coeff(-1)*zp if n==1 else coeff(2)*zm+coeff(-2)*zp+(coeff(1)**2-coeff(2))*p1**2
assert abs(r)<mp.mpf('1e-120'),(key,'Wronskian',str(r))
residual=mp.nstr(r,20)
special=None
if (N,n)==(6,2):
# Watkins Example 2.2, no hypergeometric helper.
g=(mp.sqrt(2)-1)*mp.pi**5/(32*mp.gamma(mp.mpf(3)/4)**4*mp.gamma(mp.mpf(7)/8)**2*mp.gamma(mp.mpf(5)/8)**2)
special=str(agreeing_digits(actual,g));assert mp.mpf(special)>threshold
if (N,n)==(3,2) and not stored_mode:
# Watkins Example 2.1 fixes an individual sector and its scale.
g=8*(mp.sqrt(5)-1)*mp.pi**4/(mp.mpf(5)**(mp.mpf(17)/5)*mp.gamma(mp.mpf(4)/5)**4*mp.gamma(mp.mpf(3)/5)**2)
assert agreeing_digits(row['neumann200'],g)>120
digits=18 if library else 90
values[key]=mp.nstr(mp.mpf(actual),digits)
expected[key]=mp.nstr(mp.mpf(ref),digits)
record={'key':key,'reference':ref,'full_reference':full,'agreement_digits':str(d),'stored':stored_mode,'library':library,'wronskian_residual':residual,'sextic_closed_form_digits':special}
results[key]=record
print('CHECK',json.dumps(record),flush=True)
assert not exactness(values)
assert not agrees_with(values,lambda key:expected[key])
print('SUMMARY',json.dumps({'entries':len(results),'minimum_agreement':str(min(mp.mpf(x['agreement_digits']) for x in results.values())),'stored':stored_mode,'library':library}),flush=True)
return results
if __name__=='__main__':
run_checks('--stored' in sys.argv,'--library' in sys.argv)