back to table · edit · history · where entries came from · files · download
7002 bytes, as of the version from 2026-10-29 00:57 (current). Recorded here, not run.
"""Independent numerical checks for numberdb.org/T595.
Run: sage -python checks.py > numerical-evidence.json
Keep generate.py and check.py (the campaign's reusable checks) beside this file.
This script uses mpmath at 180 decimal digits to analytically continue the
spectral sum, rather than the generator's finite Stirling-number formulas.
It is a corroboration, not the basis of the generator's proven error claim.
The explicit truncation estimate is derived in derivation.md.
"""
import json
import math
import platform
import time
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
import mpmath as mp
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.env import SAGE_VERSION
from generate import SphereDeterminants, logdet, stirling_row, even_correction
from check import exactness, measure, names_its_rings
def relative_hurwitz(p, q):
# mpmath 1.4.1's Hurwitz Euler-Maclaurin stopping test is absolute.
# For small positive values this must be converted to a relative target
# before large powers of a amplify the absolute error.
extra = max(0, int(mp.ceil(p*mp.log10(q)))) + 20
with mp.workdps(mp.mp.dps+extra):
result = mp.zeta(p, q)
return +result
def spectral_continuation(d, terms=120):
"""Direct binomial continuation of multiplicity-weighted spectral sum."""
if d == 1:
return -4*mp.diff(mp.zeta, 0), mp.mpf(0)
a = Fraction(d-1, 2)
polynomial = [Fraction(0), Fraction(2, math.factorial(d-1))]
for j in range(1, d-1):
product = [Fraction(0)]*(len(polynomial)+1)
for k, c in enumerate(polynomial):
product[k] += c*(j-a)
product[k+1] += c
polynomial = product
convert = lambda v: mp.mpf(v.numerator)/v.denominator
a = convert(a)
N = max(30, 2*d)
q = N+a
finite = mp.fsum((math.comb(l+d,d)-math.comb(l+d-2,d))
*mp.log(l*(l+d-1)) for l in range(1,N))
derivative = mp.mpf(0)
for j, coefficient in enumerate(polynomial):
if coefficient == 0:
continue
c = convert(coefficient)
part = 2*mp.zeta(-j,q,derivative=1)
for k in range(1, terms+1):
if 2*k-j == 1:
value = mp.harmonic(k-1)/2-mp.digamma(q)
else:
value = (relative_hurwitz(2*k-j,q) if 2*k-j > 1
else mp.zeta(2*k-j,q))
part += a**(2*k)*value/k
derivative += c*part
r = (a/q)**2
A = mp.fsum(abs(convert(c))*q**j for j,c in enumerate(polynomial))
bound = A*r**(terms+1)/((terms+1)*(1-r))*(1+q/(2*terms+2-d))
return finite-derivative,bound
def run_checks():
mp.mp.dps = 180
# Nonzero known controls test the continuation before comparisons.
controls = {}
expected = {1:2*mp.log(2*mp.pi),
2:mp.mpf(1)/6+4*mp.log(mp.glaisher),
3:mp.log(mp.pi)+mp.zeta(3)/(2*mp.pi**2)}
independent = {}
for d in (1,2,3):
actual,bound = spectral_continuation(d)
assert abs(actual-expected[d]) < mp.mpf('1e-140')
assert abs(actual) > mp.mpf('0.1')
independent[d] = (actual,bound)
controls[str(d)] = mp.nstr(abs(actual-expected[d]),8)
# The reusable exactness checker also sees every rational coefficient.
coeff = {str(d):list(stirling_row(d)) for d in range(21)}
coeff.update({'correction-'+str(n):even_correction(n) for n in range(1,11)})
assert not exactness(coeff)
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
g=SphereDeterminants()
values={','.join(p.values()):g.value(p,100) for p in g.enumerate()}
assert not exactness(values)
assert len(values)==40
rows=[]
C=ComplexBallField(numberdb.bits(200))
# The functional equation is an independent check of every negative-even
# zeta derivative used by the odd-dimensional finite formula.
for k in range(1,10):
lhs=C(-2*k).zetaderiv(1).real()
rhs=((-1)**k)*math.factorial(2*k)*C(2*k+1).zeta().real()/(2*(2*C.pi().real())**(2*k))
assert lhs.is_finite() and rhs.is_finite()
assert lhs.accuracy()>500 and rhs.accuracy()>500 and lhs.overlaps(rhs)
minimum_accuracy=10000
minimum_agreement=10000
table2=['39.47842','3.195311','3.338851','1.736943','1.762919',
'1.290018','1.222521','1.050409','0.946733','0.896183',
'0.778048','0.786904','0.663546','0.704655','0.580375',
'0.640108','0.517020','0.587870','0.467028','0.544581']
for d in range(1,21):
start=time.monotonic()
L150,L200=logdet(d,150),logdet(d,200)
assert L150.is_finite() and L200.is_finite() and L150.overlaps(L200)
for value in (L150,L150.exp(),L200,L200.exp()):
assert value.is_finite() and value.accuracy()>=numberdb.bits(135)
accuracy=min(L150.accuracy(),L150.exp().accuracy())
minimum_accuracy=min(minimum_accuracy,int(accuracy))
independent_L,bound=independent.get(d) or spectral_continuation(d)
Lmid=mp.mpf(str(L200.mid()))
difference=abs(Lmid-independent_L)
assert bound<mp.mpf('1e-140')
assert difference < mp.mpf('1e-130')
agreement=int(mp.floor(-mp.log10(max(difference,mp.mpf('1e-180')))))
minimum_agreement=min(minimum_agreement,agreement)
D=mp.exp(independent_L)
# Source Table 2 is rounded, so use one unit of its last place.
source=mp.mpf(table2[d-1]); ulp=mp.mpf(10)**(-len(table2[d-1].split('.')[1]))
assert abs(D-source)<ulp
rows.append({'d':d,'logdet_150_ball':str(L150),'det_150_ball':str(L150.exp()),
'accuracy_bits':int(accuracy),'independent_logdet':mp.nstr(independent_L,165),
'series_tail_bound':mp.nstr(bound,12),
'independent_absolute_difference':mp.nstr(difference,12),
'independent_agreeing_decimal_places':agreement,
'source_table_2_det':table2[d-1],'seconds':round(time.monotonic()-start,3)})
print('checked dimension',d,'accuracy bits',accuracy,'agreement places',agreement,flush=True,file=__import__('sys').stderr)
return {'table':'T595','skill_manifest':'agents/runs/20261029T002659Z-build-skills/manifest.json',
'skill_sha256':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
'versions':{'sage':SAGE_VERSION,'python':platform.python_version(),'mpmath':mp.__version__},
'controls_absolute_errors':controls,'dimensions_checked':20,'entries_checked':40,
'working_decimal_digits':[150,200],'independent_working_decimal_digits':180,
'minimum_accuracy_bits_at_150_digits':minimum_accuracy,
'minimum_independent_agreeing_decimal_places':minimum_agreement,
'measurement':measure(values),'rows':rows}
if __name__=='__main__':
print(json.dumps(run_checks(),indent=2))