back to table · edit · history · where entries came from · files · download
8045 bytes, as of the version from 2026-11-01 11:37 (current). Recorded here, not run.
"""Independent checks of T647, Broadhurst–Mellit Bessel moment determinants.
Run: sage -python checks.py > numerical-evidence.json
Place generate.py and T521.json beside this file. The T521 snapshot is a
read-only corpus comparison, not an input to the generator. Optional
T647-stored.json enables the same identities on the stored decimal intervals.
The generator uses product formulas in Arb. This file independently builds
both sequences by Zhou's coupled recursions (4.1), (4.2) in mpmath at 180
working digits; it never calls the generator to compute expected values.
Direct quadrature for every n=3..22 forms the defining matrices, with cached
Arb complex Bessel evaluations and mpmath tanh-sinh integration. Runs at
70 and 100 working digits are heuristic checks, not certification. Nonzero
controls precede them. The error certificate for the stored 100 digits is
provided by the generator's product-formula balls, not quadrature agreement.
Skill manifest: ../20261101T092014Z-build-skills/manifest.json, both hashes
a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278.
"""
import json
import time
from pathlib import Path
from functools import lru_cache
from itertools import permutations
from fractions import Fraction
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
from sage.env import SAGE_VERSION
import mpmath as mp
from generate import BroadhurstMellit, determinant
HERE = Path(__file__).resolve().parent
SKILL_HASH = 'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278'
def frac(text):
v = Fraction(text)
return QQ(v.numerator)/QQ(v.denominator)
def recorded_ball(text, R):
# Plain stored decimal means +/- one unit in its last written place.
from decimal import Decimal
d = Decimal(text)
e = d.as_tuple().exponent
u = QQ(10)**e
v = frac(text)
return R(v).add_error(R(u))
def mp_ball(x):
return mp.mpf(str(x.mid()))
def relative(a, b):
return abs((a-b)/b)
def recursive_values():
mp.mp.dps = 180
M = {1: mp.pi/(3*mp.sqrt(3))}
N = {1: mp.pi**2/16}
for k in range(2, 11):
a = mp.mpf(k)*mp.gamma(mp.mpf(k)/2)**2/(2*(2*k+1))
for j in range(1, k+1):
q = mp.mpf((2*j)**2)/((2*j)**2-1)
a *= q**(mp.mpf(k)-mp.mpf(1)/2)
M[k] = a*N[k-1]**2/M[k-1]
b = mp.mpf(2*k+1)/((k+1)*mp.factorial(k-1))
for j in range(2, k+2):
q = mp.mpf((2*j-1)**2)/((2*j-1)**2-1)
b *= q**k
N[k] = b*M[k]**2/N[k-1]
return {n: (M if n % 2 else N)[(n-1)//2] for n in range(3,23)}
def leibniz(a):
m = len(a)
total = mp.mpf(0)
for p in permutations(range(m)):
inversions = sum(p[i]>p[j] for i in range(m) for j in range(i+1,m))
term = mp.mpf((-1)**inversions)
for i in range(m):
term *= a[i][p[i]]
total += term
return total
def quadrature_run(dps, expected):
mp.mp.dps = dps
C = ComplexBallField(numberdb.bits(dps+15))
stop = 3*(dps+40)
factor_precision = [0]
@lru_cache(maxsize=40000)
def factors(t):
argument = frac(mp.nstr(t, dps+15))
precision = C.precision()
while True:
z = ComplexBallField(precision)(argument)
i = z.bessel_I(0).real()
k = z.bessel_K(0).real()
if all(v.is_finite() and v > 0 and
abs(v.mid()) > v.rad()*10**(dps+12) for v in (i,k)):
factor_precision[0] = max(factor_precision[0], int(precision))
return mp_ball(i), mp_ball(k)
precision *= 2
if precision > 10000:
raise ArithmeticError('Bessel factors did not attain target accuracy')
def moment(a,b,r):
def f(t):
if t == 0:
return mp.mpf(0)
i,k = factors(t)
return i**a*k**b*t**r
return mp.quad(f,[0,1,5,20,stop],method='tanh-sinh')
controls=[]
for a,b,r,known in [(0,1,1,mp.mpf(1)),(0,2,1,mp.mpf(1)/2)]:
q=moment(a,b,r);err=relative(q,known)
assert err < mp.mpf(10)**(-dps+10), ('control',a,b,r,mp.nstr(q,dps),mp.nstr(err,20))
controls.append({'a':a,'b':b,'r':r,'value':mp.nstr(q,dps),'relative_error':mp.nstr(err,8)})
rows=[]
for n in range(3,23):
m=(n-1)//2
a=[[moment(i,n-i,2*j-1) for j in range(1,m+1)] for i in range(1,m+1)]
value=mp.det(mp.matrix(a))
err=relative(value,expected[n])
assert err < mp.mpf('1e-40'), (n,dps,mp.nstr(err))
if m<=3:
assert relative(leibniz(a),value) < mp.mpf(10)**(-dps+10)
# Ratio of product of row 1-norms to determinant is a conservative
# size/cancellation diagnostic, not an integration error estimate.
diagnostic=mp.log10(mp.fprod(sum(abs(x) for x in row) for row in a)/abs(value))
rows.append({'n':n,'determinant':mp.nstr(value,dps-10),'relative_error':mp.nstr(err,8),
'log10_row_norm_product_over_det':mp.nstr(diagnostic,8)})
print('quadrature',dps,n,'relative error',mp.nstr(err,5),file=__import__('sys').stderr,flush=True)
return {'working_digits':dps,'cutoff':stop,'controls':controls,'rows':rows,'cached_nodes':factors.cache_info().currsize,'max_bessel_working_bits':factor_precision[0]}
def main():
started=time.time()
expected=recursive_values()
R=RealBallField(800)
g=BroadhurstMellit()
computed={int(p['n']):g.value(p,100) for p in g.enumerate()}
rows=[]
for n,v in computed.items():
high=determinant(n,180)
assert v.is_finite() and v.overlaps(high) and not v.contains_zero()
err=relative(mp_ball(high),expected[n])
assert err < mp.mpf('1e-165'),(n,err)
text=to_text(v,digits=100)
assert R(high) in recorded_ball(text,R)
radius=mp.mpf(str(v.rad()))
accuracy=-mp.log10(radius/abs(mp_ball(v)))
rows.append({'n':n,'number':text,'relative_radius_decimal_digits':mp.nstr(accuracy,10),
'recursion_relative_error':mp.nstr(err,10)})
# An intentionally wrong value must fail the interval check.
control=recorded_ball(rows[0]['number'],R)
assert not control.overlaps(R(computed[3])*R(1+QQ(1)/QQ(10)**90))
corpus=json.loads((HERE/'T521.json').read_text())['Numbers']
corpus_checks=[]
for n in range(3,7):
m=(n-1)//2
a=[[mp.mpf(corpus[f'{i},{n-i}'][str(2*j-1)]) for j in range(1,m+1)] for i in range(1,m+1)]
q=leibniz(a)
err=relative(q,expected[n]); assert err<mp.mpf('1e-97')
corpus_checks.append({'n':n,'relative_error':mp.nstr(err,8)})
quadratures=[quadrature_run(p,expected) for p in (70,100)]
stored_checks=[]
path=HERE/'T647-stored.json'
if path.exists():
mp.mp.dps=180
doc=json.loads(path.read_text()); numbers=doc['Numbers']
assert list(map(int,numbers))==list(range(3,23))
for n in range(3,23):
text=numbers[str(n)]
if isinstance(text,dict):text=text['number']
b=recorded_ball(text,R)
assert R(determinant(n,180)) in b
# mpmath recurrence approximation is checked inside the written
# interval with much more precision than the stored last place.
assert R(frac(mp.nstr(expected[n],170))) in b
stored_checks.append(n)
print('JSON_RESULT=' + json.dumps({'table':'T647','skill_manifest':'../20261101T092014Z-build-skills/manifest.json',
'skill_sha256':{'local':SKILL_HASH,'live':SKILL_HASH},
'versions':{'sage':SAGE_VERSION,'numberdb':numberdb.__version__,'mpmath':mp.__version__},
'product_ball_checks':rows,'corpus_determinants':corpus_checks,'quadratures':quadratures,
'stored_interval_checks':stored_checks,'negative_control_passed':True,'elapsed_seconds':time.time()-started},indent=2))
if __name__=='__main__':
main()