back to table · edit · history · where entries came from · files · download
4886 bytes, as of the version from 2026-10-16 07:34 (current). Recorded here, not run.
"""Independent repair checks for Kolmogorov quantiles -- numberdb.org/T456.
Run with SageMath and the downloaded table document:
sage -pip install numberdb
sage -python repair_checks.py < T456-before.json
Use 600-bit MPFI and Jacobi's product (DLMF 20.5.4), independent of the
attached generator's alternating-series inversion. With q=exp(-2*x*x),
P_N=product((1-q**(2*k))*(1-q**(2*k-1))**2, k=1..N).
For the omitted factors, sum their deficits is
S=(q**(2*N+2)+2*q**(2*N+1))/(1-q*q).
For 0<q<1 and S<1, P_N*(1-S) <= theta4(0,q) <= P_N.
Every arithmetic operation is an outward-rounded interval operation. We
require finite enclosures of diameter <1e-170, not merely overlapping balls.
The stored decimal endpoints are parsed exactly and bracket F(x)=p.
The identity added to Similar tables follows term by term from T426:
theta4(0,q)=1+2*sum((-1)**j*q**(j*j)); substitute q=exp(-2*x*x).
For x>0 this nome lies in (0,1), and the series converges absolutely.
The relation to the scaled statistic is the limiting theorem stated by SciPy;
the checks confirm the alpha=1-p tail convention for every stored row.
"""
import ast
import json
import sys
from decimal import Decimal
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.real_arb import RealBallField
import mpmath as mp
from scipy.stats import kstwobign
I = RealIntervalField(600)
R = RealBallField(600)
N = 600
def rational_decimal(s):
d = Decimal(s)
f = Fraction(d)
return QQ(f.numerator)/f.denominator, QQ(10)**d.as_tuple().exponent
def narrow(v):
assert not v.is_NaN() and not v.lower().is_infinity() and not v.upper().is_infinity()
assert I(v.absolute_diameter()) < I(10)**-170
return v
def product_cdf(x):
x = I(x)
assert x > 0
q = (-2*x*x).exp()
assert q > 0 and q < 1
product = I(1)
for k in range(1, N+1):
product *= (1-q**(2*k))*(1-q**(2*k-1))**2
deficit = (q**(2*N+2)+2*q**(2*N+1))/(1-q*q)
assert deficit > 0 and deficit < 1
return narrow(product*I(1-deficit.upper(), 1))
def run_program(source):
tree = ast.parse(source)
assert isinstance(tree.body[-1], ast.Expr)
ns = {}
exec(compile(ast.Module(body=tree.body[:-1], type_ignores=[]), '<table-program>', 'exec'), ns)
return eval(compile(ast.Expression(tree.body[-1].value), '<table-program>', 'eval'), ns)
doc = json.load(sys.stdin)
control = product_cdf(1)
assert control > I(QQ(72)/100) and control < I(QQ(74)/100)
# A special-value control independent of the product/series representations.
special = product_cdf((I.pi()/2).sqrt())
gamma = narrow(I(R(QQ(1)/4).gamma()/(2*R.pi())**(QQ(3)/4)))
assert special.overlaps(gamma)
assert kstwobign.cdf(0) == 0 and kstwobign.sf(0) == 1
mp.mp.dps = 150
rows = []
previous = None
for p, text in doc['Numbers'].items():
midpoint, ulp = rational_decimal(text)
left, right = product_cdf(midpoint-ulp), product_cdf(midpoint+ulp)
target = I(QQ(p))
assert left.upper() < target.lower() and right.lower() > target.upper(), p
if previous is not None:
assert midpoint-ulp > previous
previous = midpoint+ulp
assert product_cdf(midpoint+3*ulp).lower() > target.upper()
assert product_cdf(midpoint-3*ulp).upper() < target.lower()
alpha = I(1)-target
assert (1-right).upper() < alpha.lower() and (1-left).lower() > alpha.upper()
x = mp.mpf(str(midpoint.numerator()))/mp.mpf(str(midpoint.denominator()))
theta = mp.jtheta(4, 0, mp.exp(-2*x*x))
direct = 1-2*mp.fsum((-1)**(j-1)*mp.exp(-2*j*j*x*x) for j in range(1,81))
discrepancy = abs(theta-direct)
assert discrepancy < mp.mpf('1e-140')
assert abs(float(kstwobign.ppf(float(QQ(p))))-float(midpoint)) < 2e-12
rows.append({'p':p, 'product_bracket':True, 'shift_controls_rejected':True,
'upper_tail_without_halving':True,
'theta_series_discrepancy':mp.nstr(discrepancy,8),
'cdf_width_bound':'1e-170'})
programs = {}
midpoint, ulp = rational_decimal(doc['Numbers']['19/20'])
for name, program in doc['Programs'].items():
value = run_program(program['code'])
if name == 'program-mpmath':
reference = mp.mpf(str(midpoint.numerator()))/mp.mpf(str(midpoint.denominator()))
scale = mp.mpf(str(ulp.numerator()))/mp.mpf(str(ulp.denominator()))
error = abs(value-reference)/scale
assert error < 1
programs[name] = {'error_in_last_place':mp.nstr(error,12)}
else:
error = abs(float(value)-float(midpoint))
assert error < 2e-12
programs[name] = {'absolute_error':error}
print(json.dumps({'table':'T456','ok':True,'checked':len(rows),'bits':600,
'product_terms':N,'known_controls':True,
'entries':rows,'programs':programs},indent=2))