back to table · edit · history · where entries came from · files · download
8475 bytes, as of the version from 2026-10-24 21:16 (current). Recorded here, not run.
"""Independent integral checks for numberdb.org/T467.
sage -python checks.py # no key or database required
sage -python checks.py --stored # also test database values
Requires mpmath, SageMath/passagemath, numberdb, and the attached generate.py
and check.py. A draft read requires NUMBERDB_API_KEY (or a key on stdin with
NUMBERDB_KEY_FROM_STDIN=1). NUMBERDB_CHECK_STORED=1 also enables read-back.
The independent computations use mpmath at 130 and 170 decimal digits and
must agree with the 100-digit intervals. Their error is not certified; the
proven claim comes from Arb evaluation of published closed forms.
Box method: the absolute difference of two independent U[0,L] coordinates
has density 2(L-x)/L^2. Its Gaussian transform is
F(L*t)=sqrt(pi)*erf(L*t)/(L*t)-(1-exp(-(L*t)^2))/(L*t)^2.
The identity sqrt(s)=integral_0^infinity(1-exp(-s*t^2))/t^2 dt/sqrt(pi)
reduces each box mean to a one-dimensional integral of 1-product F(L*t).
For t<1 we form 1-F by its convergent series to avoid cancellation. For
t>1 substitute t=1/u, so no infinite tail is dropped.
Equilateral method: the difference body is a regular hexagon. On a ray
z=t*w to an edge point w, the intersection area is A*(1-t)^2. Integrating
t^2*(1-t)^2 gives 1/30. Six edge determinants divided by A give 12, hence
the mean is (2/5)*integral_0^1 sqrt(1-u+u^2) du. The zeroth moment is 1.
The disk uses the overlap density from Tu/Fischbach (12), with r=2 sin(t).
The solid ball uses their polynomial density (15), integrated exactly.
Circle/sphere use their uniform angle / height coordinate distributions.
Every numerical method first checks a known answer (normalization and a
second moment where available). No Monte Carlo check supplies digits here.
"""
import json
import os
import sys
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
import mpmath as mp
from generate import MeanDistance, REGIONS, closed_form
from check import exactness, measure, names_its_rings, agrees_with
from numberdb._write import to_text
def close(x, y, digits=115):
assert mp.isfinite(x) and mp.isfinite(y)
assert abs(x-y) < mp.mpf(10)**(-digits), (mp.nstr(x,30),mp.nstr(y,30))
def box_mean(lengths):
def deficit(z):
# 1-F(z), summing exact moment coefficients in arbitrary precision.
term=z*z/6
out=term
k=1
while term:
term *= -z*z*(2*k+1)*(2*k+2)/((k+1)*(2*k+3)*(2*k+4))
out += term
k += 1
if abs(term) < mp.eps*abs(out):
return out
assert k < 3000
return out
def low(t):
if not t:
return sum(mp.mpf(a)**2 for a in lengths)/6
d=mp.mpf(0)
for a in lengths:
e=deficit(a*t)
d=d+e-d*e
return d/(t*t)
def high(u):
if not u:return mp.mpf(1)
product=mp.mpf(1)
for a in lengths:
z=mp.mpf(a)/u
product *= mp.sqrt(mp.pi)*mp.erf(z)/z+mp.expm1(-z*z)/(z*z)
return 1-product
return (mp.quad(low,[0,1])+mp.quad(high,[0,1]))/mp.sqrt(mp.pi)
def independent(dps):
mp.mp.dps=dps
# Controls precede the computations that rely on them.
for length in (1,2):
L=mp.mpf(length)
density=lambda x:2*(L-x)/(L*L)
close(mp.quad(density,[0,L]),1)
close(mp.quad(lambda x:x*x*density(x),[0,L]),L*L/6)
close(box_mean((1,)),mp.mpf(1)/3)
close(box_mean((2,)),mp.mpf(2)/3)
out={'unit-interval':mp.mpf(1)/3}
for region,lengths in [('unit-square',(1,1)),('rectangle-1x2',(1,2)),('unit-cube',(1,1,1))]:
out[region]=box_mean(lengths)
close(12*mp.quad(lambda t:t*(1-t)**2,[0,1]),1)
out['equilateral-triangle']=mp.mpf(2)/5*mp.quad(lambda u:mp.sqrt(1-u+u*u),[0,1])
# Disk density times dr/dt under r=2 sin(t), 0<=t<=pi/2.
def disk_weight(t):
return 16/mp.pi*mp.sin(t)*mp.cos(t)*(mp.pi/2-t-mp.sin(t)*mp.cos(t))
close(mp.quad(disk_weight,[0,mp.pi/2]),1)
close(mp.quad(lambda t:(2*mp.sin(t))**2*disk_weight(t),[0,mp.pi/2]),1)
out['unit-disk']=mp.quad(lambda t:2*mp.sin(t)*disk_weight(t),[0,mp.pi/2])
# Exact antiderivatives of p(r)=3r^2-9r^3/4+3r^5/16, 0<=r<=2.
def ball_moment(m):
return sum(c*QQ(2)**(k+m+1)/(k+m+1)
for k,c in [(2,QQ(3)),(3,-QQ(9)/4),(5,QQ(3)/16)])
assert ball_moment(0)==1 and ball_moment(2)==QQ(6)/5
assert not agrees_with({'ball-3':ball_moment(1)},lambda key:QQ(36)/35)
out['unit-ball-3']=mp.mpf(str(ball_moment(1).numerator()))/int(ball_moment(1).denominator())
close(mp.quad(lambda t:mp.mpf(1)/mp.pi,[0,mp.pi]),1)
close(mp.quad(lambda t:(2*mp.sin(t/2))**2/mp.pi,[0,mp.pi]),2)
out['unit-circle']=mp.quad(lambda t:2*mp.sin(t/2)/mp.pi,[0,mp.pi])
# On the sphere, z is uniform on [-1,1] and the chord length is sqrt(2-2z).
close(mp.quad(lambda z:mp.mpf(1)/2,[-1,1]),1)
close(mp.quad(lambda z:(2-2*z)/2,[-1,1]),2)
out['unit-sphere-2']=mp.quad(lambda z:mp.sqrt(2-2*z)/2,[-1,1])
# Scaling checks on every method / region, with scale 2.
close(box_mean((2,2)),2*out['unit-square'])
close(box_mean((2,4)),2*out['rectangle-1x2'])
close(box_mean((2,2,2)),2*out['unit-cube'])
close(mp.mpf(2)/5*mp.quad(lambda u:mp.sqrt(4-4*u+4*u*u),[0,1]),2*out['equilateral-triangle'])
close(mp.quad(lambda t:4*mp.sin(t)*disk_weight(t),[0,mp.pi/2]),2*out['unit-disk'])
assert 2*ball_moment(1)==QQ(72)/35
close(mp.quad(lambda t:4*mp.sin(t/2)/mp.pi,[0,mp.pi]),2*out['unit-circle'])
close(mp.quad(lambda z:mp.sqrt(8-8*z)/2,[-1,1]),2*out['unit-sphere-2'])
print('controls and scaling passed at',dps,'digits',flush=True)
print('INTEGRALS_'+str(dps),json.dumps({k:mp.nstr(out[k],125) for k in REGIONS}),flush=True)
return out
def mp_number(text):
if '/' in text:
a,b=text.split('/')
return mp.mpf(a)/mp.mpf(b)
return mp.mpf(text)
def interval_of(text):
text=str(text).strip()
if '/' in text or ('.' not in text and 'e' not in text.lower()):
v=mp_number(text)
return v,v
if '+/-' in text:
a,b=text.split('+/-');v=mp.mpf(a);r=mp.mpf(b);return v-r,v+r
if text.startswith('['):
a,b=text.strip('[]').split(',');return mp.mpf(a),mp.mpf(b)
mantissa,*exponent=text.lower().split('e')
decimals=len(mantissa.split('.')[1]) if '.' in mantissa else 0
place=(int(exponent[0]) if exponent else 0)-decimals
v=mp.mpf(text);r=mp.mpf(10)**place
return v-r,v+r
def check_text(text, v):
lo,hi=interval_of(text)
assert mp.isfinite(lo) and mp.isfinite(hi)
if lo==hi:close(v,lo)
else:
assert hi-lo < mp.mpf('3e-99')
assert lo < v < hi,(text,mp.nstr(v,110))
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN'):
numberdb.configure(api_key=sys.stdin.read().strip())
first=independent(130)
second=independent(170)
g=MeanDistance()
values={k:g.value({'region':k},100) for k in REGIONS}
assert not exactness(values)
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
serialized={}
radii={}
for k in REGIONS:
close(first[k],second[k],120)
v=values[k]['number']
text=to_text(v,digits=100)
serialized[k]=text
check_text(text,second[k])
if hasattr(v,'rad'):
assert v.is_finite() and v.rad() < QQ(1)/10**120
radii[k]=str(RealField(40)(v.rad()))
# Higher-precision Arb calculation must refine consistently.
higher=closed_form(k,150)
assert higher.is_finite() and v.overlaps(higher)
print(k,text,flush=True)
print('SERIALIZED_MEASURE',json.dumps(measure(serialized)),flush=True)
print('BALL_RADII',json.dumps(radii),flush=True)
print('INDEPENDENT_VALUES',json.dumps({k:mp.nstr(second[k],135) for k in REGIONS}),flush=True)
if '--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1':
document=numberdb.table('T467')
from numberdb._generate import _stored_entries
entries=_stored_entries(document)
assert len(entries)==len(REGIONS)
for k,text in entries.items():
check_text(text,second[k])
report=g.verify(sample=None)
assert report.ok,report
print('STORED_CHECK',report,flush=True)
print('ALL_CHECKS_PASSED',flush=True)
if __name__=='__main__':main()