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