back to table · edit · history · where entries came from · files · download
6878 bytes, as of the version from 2026-10-16 07:34 (current). Recorded here, not run.
"""Independent checks for Kolmogorov quantiles -- numberdb.org/T456.
sage -python checks.py # recompute and check before writing
sage -python checks.py --stored # check values read back from NumberDB
The independent CDF uses MPFI interval arithmetic, the positive theta series,
and a geometric tail bound. It shares neither the summation nor the rounding
backend with the Arb alternating-series generator. After N terms, put
a = pi^2/(8*x^2), b = exp(-a*(2*N+1)^2), rho = exp(-8*a*(N+1)).
Successive omitted terms have ratio at most rho, so the omitted sum is at
most b/(1-rho). Each decimal in NumberDB means its centre plus/minus one unit
in the last place. The check proves F(lower)<p<F(upper) for those endpoints.
Controls: F(0)=0; at x=sqrt(pi/2), F(x)=Gamma(1/4)/(2*pi)^(3/4).
The latter follows from theta2(i)=theta4(i), theta3^4=theta2^4+theta4^4,
and K(1/sqrt(2))=Gamma(1/4)^2/(4*sqrt(pi)). It checks the scaling and
prefactor of both CDF series against an independent special value.
Every interval is checked for finiteness and a meaningful width.
SciPy comparisons are deliberately limited to 12 decimal places.
"""
import json
import os
import sys
from decimal import Decimal
from fractions import Fraction
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.real_arb import RealBallField
from scipy.stats import kstwobign
import scipy
import mpmath as mp
from generate import KolmogorovQuantiles, PROBABILITIES, cdf
from check import exactness, measure, names_its_rings
CHECK_BITS = 560
I = RealIntervalField(CHECK_BITS)
def transformed_cdf(x):
"""Positive theta series with a rigorous geometric remainder in MPFI."""
x = I(x)
if x <= 0:
return I(0)
assert x > 0
a = I.pi()**2/(8*x*x)
total = I(0)
prefactor = (2*I.pi()).sqrt()/x
for k in range(1, 10000):
total += (-a*(2*k-1)**2).exp()
first_omitted = (-a*(2*k+1)**2).exp()
ratio = (-8*a*(k+1)).exp()
tail = first_omitted/(1-ratio)
if I(prefactor.upper()*tail.upper()) < I(10)**(-145):
ans = prefactor*(total+I(0,tail.upper()))
assert not ans.is_NaN() and ans.lower() > -1 and ans.upper() < 2
assert I(ans.absolute_diameter()) < I(10)**(-140)
return ans
raise ArithmeticError('independent theta series did not converge')
def decimal_interval(text):
"""Parse a stored decimal without asking NumberDB's parser."""
d = Decimal(text)
f = Fraction(d)
mid = QQ(f.numerator)/QQ(f.denominator)
exponent = d.as_tuple().exponent
ulp = QQ(10)**exponent
return mid-ulp, mid+ulp, mid, ulp
def control():
assert transformed_cdf(0) == 0
R = RealBallField(CHECK_BITS)
x = (I.pi()/2).sqrt()
known_ball = R(QQ(1)/4).gamma()/(2*R.pi())**(QQ(3)/4)
known = I(known_ball)
theta = transformed_cdf(x)
direct = I(cdf((R.pi()/2).sqrt(),R))
for value in [known, theta, direct]:
assert not value.is_NaN()
assert I(value.absolute_diameter()) < I(10)**(-135)
assert theta.overlaps(known) and direct.overlaps(known)
# The reference API's survival and lower-tail conventions at zero.
assert kstwobign.cdf(0) == 0 and kstwobign.sf(0) == 1
return {'gamma_control_value': to_text(known_ball, digits=100),
'theta_gamma_overlap': True, 'alternating_gamma_overlap': True,
'width_bound': '1e-135', 'F(0)': 0}
def mpmath_quantile(p):
"""Reader-facing snippet, independent special-function implementation."""
mp.mp.dps = 130
a,b = str(p).split('/')
probability = mp.mpf(a)/mp.mpf(b)
return mp.findroot(lambda x: mp.jtheta(4,0,mp.exp(-2*x*x))-probability,
(mp.mpf('0.25'),mp.mpf(4)),solver='bisect',
maxsteps=600, tol=mp.mpf('1e-120'))
def run(stored=False):
evidence = {'table':'T456', 'source':'database' if stored else 'generator',
'controls':control(), 'check_bits':CHECK_BITS,
'scipy_version':scipy.__version__, 'entries':[]}
g=KolmogorovQuantiles()
if stored:
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
doc=numberdb.table(g.table)
records=doc['Numbers']
assert isinstance(records,dict), type(records)
texts={str(p):v['number'] if isinstance(v,dict) else v
for p,v in records.items()}
assert tuple(texts) == PROBABILITIES
else:
values={p:g.value({'p':p},100) for p in PROBABILITIES}
assert exactness(values)==[]
assert names_its_rings(__file__.replace('checks.py','generate.py'))==[]
texts={p:to_text(v,digits=100) for p,v in values.items()}
evidence['raw_size']=measure(values)
evidence['max_ball_radius']=str(max(v.rad() for v in values.values()))
evidence['minimum_retained_decimal_digits']=min(
len(to_text(v,digits=130).replace('.','').lstrip('0'))
for v in values.values())
previous=None
for p,text in texts.items():
low,high,mid,ulp=decimal_interval(text)
target=I(QQ(p))
left,right=transformed_cdf(low),transformed_cdf(high)
assert QQ(left.upper()) < QQ(target.lower()), ('lower sign',p)
assert QQ(right.lower()) > QQ(target.upper()), ('upper sign',p)
if previous is not None:
assert low>previous, ('monotonicity',p)
previous=high
# Demonstrate that a corrupted written interval is rejected.
shifted=transformed_cdf(low+100*ulp)
assert QQ(shifted.lower())>QQ(target.upper())
fref=float(kstwobign.ppf(float(QQ(p))))
err=abs(fref-float(mid))
assert err<2e-12, ('scipy',p,err)
independent=mpmath_quantile(p)
exact_mp=mp.mpf(str(mid.numerator()))/mp.mpf(str(mid.denominator()))
err_mp=abs(independent-exact_mp)
mp_ulp=mp.mpf(str(ulp.numerator()))/mp.mpf(str(ulp.denominator()))
assert err_mp<mp_ulp, ('mpmath',p)
# Upper-tail conversion on every row, with no division by two.
alpha=I(1)-target
assert QQ((I(1)-right).upper())<QQ(alpha.lower())
assert QQ((I(1)-left).lower())>QQ(alpha.upper())
evidence['entries'].append({'p':p,'number':text,'theta_bracket':True,
'shifted_interval_rejected':True,'upper_tail_conversion':True,
'scipy_ppf':fref,'scipy_absolute_error':err,
'mpmath_error_in_last_place':mp.nstr(err_mp/mp_ulp,12)})
evidence['written_size']=measure(texts)
evidence['checked']=len(texts)
evidence['ok']=True
print(json.dumps(evidence,indent=2))
return evidence
if __name__=='__main__':
run(stored='--stored' in sys.argv)