"""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)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()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)