back to table · edit · history · where entries came from · files · download
7119 bytes, as of the version from 2026-10-16 07:53. Recorded here, not run.
"""Independent numerical checks for numberdb.org/T457.
Run sage -python checks.py to check generated values. For database values:
NUMBERDB_KEY_FROM_STDIN=1 sage -python checks.py --stored < key-file
The key is used only to read the private NumberDB draft.
The check uses finite elementary integration recurrences, independently of
both the generator's beta evaluation and its infinite coefficient recurrence.
"""
import json
import os
import sys
from decimal import Decimal
from time import time
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from sage.arith.misc import binomial
from scipy.stats import studentized_range
B=RealBallField(2048)
def cos_integral(n, theta):
"""Integral_0^theta cos(u)^n du by elementary integration by parts."""
s,c=theta.sin(),theta.cos()
vals=[theta,s]
for j in range(2,n+1):
vals.append(s*c**(j-1)/j+B(j-1)/j*vals[j-2])
return vals[n]
def independent_cdf(q,k,nu):
x=B(q)
if k==2:
theta=(x/B(2*nu).sqrt()).arctan()
return cos_integral(nu-1,theta)/cos_integral(nu-1,B.pi()/2)
a=x*x/(2*nu)
A=1+a
if nu%2:
m=(nu-1)//2
u=(1/(2*A.sqrt())).arcsin()
z=u.tan()
ints=[u]
for j in range(1,m+1):
ints.append(z**(2*j-1)/(2*j-1)-ints[-1])
integral=sum(B(binomial(m,j))*A**(-m)*(-a)**j*ints[j]
for j in range(m+1))
else:
m=nu//2
u=(a/(3*A)).sqrt().arctan()
# Partial fractions: I_m = I_(m-1) - a K_m, I_0=pi/6.
# Taking 1-6 I_m/pi cancels the initial constant algebraically.
return 6/B.pi()*(a*A).sqrt()*sum(
cos_integral(2*j-2,u)/A**j for j in range(1,m+1))
return 1-6/B.pi()*integral
def endpoints(text):
# The database decimal denotes +/- one unit in its last written place.
d=Decimal(text)
tup=d.as_tuple()
signed=int(''.join(str(x) for x in tup.digits))*(-1 if tup.sign else 1)
unit=QQ(10)**int(tup.exponent)
center=QQ(signed)*unit
return center-unit,center+unit,center,unit
def controls():
tests=[(independent_cdf(B(2).sqrt(),2,1),B(1)/2),
(independent_cdf(2/B(3).sqrt(),2,2),B(1)/2),
(independent_cdf((2+2*B(3).sqrt()).sqrt(),3,1),B(1)/2)]
for got,expected in tests:
assert got.is_finite() and B(got.rad())<B(10)**(-500)
assert (got-expected).contains_zero(),(got,expected)
# A nontrivial control of the angular formula against the general CDF.
got=independent_cdf(3,3,10)
scipy=studentized_range.cdf(3,3,10)
assert abs(float(got)-scipy)<1e-12
return len(tests)+1
def records_from_document(doc):
if isinstance(doc['Numbers'],list):
return [(entry['params'],entry['number']) for entry in doc['Numbers']]
records=[]
for nu,by_k in doc['Numbers'].items():
for k,by_p in by_k.items():
for p,entry in by_p.items():
text=entry['number'] if isinstance(entry,dict) else entry
records.append(({'nu':nu,'k':k,'p':p},text))
return records
def run(records):
control_count=controls()
worst_scipy=0
worst_at=None
low_margin=None
min_accuracy=None
by_params={}
closed_forms=0
started=time()
for index,(params,text) in enumerate(records,1):
nu,k,p=int(params['nu']),int(params['k']),QQ(params['p'])
low,high,center,unit=endpoints(text)
f_low=independent_cdf(low,k,nu)
f_high=independent_cdf(high,k,nu)
assert f_low.is_finite() and f_high.is_finite()
assert f_low<p and f_high>p,(params,text,f_low,f_high)
# The upper-tail interpretation is the same sign reversal at 1-p.
assert 1-f_low>1-p and 1-f_high<1-p
# A deliberately shifted last-place interval must be rejected.
shifted=independent_cdf(high+unit,k,nu)
assert shifted>p,(params,'negative control failed')
acc=min(f_low.accuracy(),f_high.accuracy())
min_accuracy=acc if min_accuracy is None else min(acc,min_accuracy)
sc=float(studentized_range.cdf(float(center),k,nu))
err=abs(sc-float(p))
if err>worst_scipy:worst_scipy,worst_at=err,params
assert err<2e-10,(params,err)
if nu==1:
if k==2: closed=B(2).sqrt()*(B.pi()*B(p)/2).tan()
else: closed=(1/(2*(B.pi()*(1-B(p))/6).sin()**2)-2).sqrt()
assert closed>low and closed<high
closed_forms+=1
elif nu==2 and k==2:
closed=2*B(p)/(1-B(p)**2).sqrt()
assert closed>low and closed<high
closed_forms+=1
by_params[nu,k,str(p)]=(low,high)
if index%56==0: print('checked',index,'elapsed',round(time()-started,2),flush=True)
# Strict monotonicity in p and the stochastic ordering in k.
from generate import DEGREES,PROBABILITIES
for nu in DEGREES:
for k in (2,3):
for p,q in zip(PROBABILITIES,PROBABILITIES[1:]):
assert by_params[nu,k,p][1]<by_params[nu,k,q][0]
for p in PROBABILITIES:
assert by_params[nu,2,p][1]<by_params[nu,3,p][0]
printed_rows=json.loads(Path(__file__).with_name('reference-values.json').read_text())
printed_rounding_discrepancies=[]
for ref in printed_rows:
low,high=by_params[ref['nu'],ref['k'],ref['p']]
rl,rh,rc,unit=endpoints(ref['printed'])
# Record rounding discrepancies; never hide them by a wider tolerance.
# One printed value is not correctly rounded: 90.03 instead of 90.02.
if not (low>rc-unit/2 and high<rc+unit/2):
printed_rounding_discrepancies.append(ref)
assert ref=={'nu':1,'k':2,'p':'99/100','printed':'90.03'}
assert low>rc-unit and high<rc+unit,ref
report={'entries':len(records),'controls':control_count,'closed_form_checks':closed_forms,
'printed_reference_checks':len(printed_rows),
'printed_rounding_discrepancies':printed_rounding_discrepancies,
'independent_endpoint_brackets':len(records),
'shifted_interval_rejections':len(records),
'upper_tail_checks':len(records),
'scipy_cdf_checks':len(records),
'worst_scipy_cdf_absolute_error':worst_scipy,
'worst_scipy_at':worst_at,'minimum_check_accuracy_bits':int(min_accuracy),
'monotonicity_in_p':True,'stochastic_order_in_k':True,
'elapsed_seconds':round(time()-started,2)}
print('CHECK_REPORT '+json.dumps(report),flush=True)
return report
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
if '--stored' in sys.argv or os.environ.get('CHECK_STORED')=='1':
records=records_from_document(numberdb.table('T457'))
else:
from generate import StudentizedRangeQuantiles
from numberdb._write import to_text
g=StudentizedRangeQuantiles()
records=[(p,to_text(g.value(p,100),100)) for p in g.enumerate()]
run(records)