back to table · edit · history · where entries came from · files · download
6672 bytes, as of the version from 2026-11-09 09:32 (current). Recorded here, not run.
"""Validation for Higher Mahler measures -- numberdb.org/T686.
Run: sage -python checks.py, with generate.py beside this file.
Add --stored with NUMBERDB_API_KEY configured to check database readback.
The radial three-step quadrature and four-step closed form and pointwise
hypergeometric density checks are independent of the generator's recurrences.
"""
import sys,os,json,math
from functools import lru_cache
from time import monotonic
import numberdb.sage as numberdb
from generate import computed,density4_series,HigherMahler
import mpmath as mp
MANIFEST='agents/runs/20261109T084515Z-build-skills/manifest.json'
SKILL_SHA256='886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b'
evidence={'skill_manifest':MANIFEST,'local_skill_sha256':SKILL_SHA256,'live_skill_sha256':SKILL_SHA256,'checks':{}}
mp.mp.dps=240
def mpnum(x):return mp.mpf(x.mid().str(truncate=False))
def compare(name,a,b,digits):
a,b=mp.mpf(a),mp.mpf(b)
err=abs(a-b)/max(1,abs(b))
assert mp.isfinite(a) and mp.isfinite(b) and err<mp.power(10,-digits),(name,mp.nstr(err,10))
evidence['checks'][name]={'relative_to_max_one_error':mp.nstr(err,12),'required_decimal_places':digits}
print('PASS',name,'error',mp.nstr(err,8),flush=True)
@lru_cache(None)
def p3(x):
if x==0 or x==1:return mp.mpf(0) # isolated quadrature endpoint values
if x<1:
d=(1+x)**3*(3-x)
complement=(1-x)**3*(x+3)/d
return 4*x/(mp.pi**2*mp.sqrt(d))*mp.elliprf(0,complement,1)
complement=(x-1)**3*(x+3)/(16*x)
return mp.sqrt(x)/mp.pi**2*mp.elliprf(0,complement,1)
def radial3(k):
def f(x):
if x==0 or x==1:return mp.mpf(0)
return p3(x)*mp.log(x)**k
return mp.quad(f,[0,mp.mpf(1)/2,1,2,3])
vals={n:[mp.mpf(s) for s in computed(n,220)] for n in (3,4)}
for n in (3,4):
for k,s in enumerate(computed(n,150)):
compare('precision_%s_%s'%(n,k),s,vals[n][k],148)
mp.mp.dps=130
# Independent controls before believing the quadrature.
compare('radial3_normalization',radial3(0),1,115)
compare('radial3_squared_radius',mp.quad(lambda x:p3(x)*x*x,[0,mp.mpf(1)/2,1,2,3]),3,115)
cl2=mp.im(mp.polylog(2,mp.exp(mp.j*mp.pi/3)))
compare('radial3_first_moment',radial3(1),cl2/mp.pi,115)
for k in range(2,9):compare('radial3_order_%d'%k,radial3(k),vals[3][k],115)
compare('three_first_moment',vals[3][1],cl2/mp.pi,120)
compare('four_first_moment',vals[4][1],7*mp.zeta(3)/(2*mp.pi**2),120)
l=mp.log(2)
closed4=(24*mp.polylog(4,mp.mpf(1)/2)-18*mp.zeta(4)+21*mp.zeta(3)*l-6*mp.zeta(2)*l*l+l**4)/mp.pi**2
compare('four_second_moment_BSWZ_70',vals[4][2],closed4,120)
for k,s in {2:'0.41929927830117445534618570174886146566170299117521',3:'0.13072798584098927059592540295887788768895327503289',4:'0.52153569858138778267996782141801173128244973155094',5:'-0.46811264825699083401802243892432823881642492433794'}.items():
compare('published_BBSW_%d'%k,vals[3][k],s,49)
# Independently evaluated hypergeometric density across every series region.
R,P,N,r,s,C,S,F=density4_series(150)
def series_density(x):
xx=R(str(x))
if x<1:return sum((r[j]-s[j]*xx.log())*xx**(2*j+1) for j in range(len(r)))
if x<=2:return C(xx-2)
if x<=3:return C(xx-2)+(xx-2).sqrt()*S(xx-2)
return (4-xx).sqrt()*F(4-xx)
for x in [mp.mpf(1)/2,mp.mpf(3)/2,mp.mpf(5)/2,mp.mpf(7)/2]:
z=(16-x*x)**3/(108*x**4)
direct=2/mp.pi**2*mp.sqrt(16-x*x)/x*mp.re(mp.hyper([mp.mpf(1)/2]*3,[mp.mpf(5)/6,mp.mpf(7)/6],z))
compare('four_density_'+str(x),mpnum(series_density(x)),direct,115)
compare('endpoint_coefficient_sign',mpnum(F[1]/F[0]),mp.mpf(3)/16,120)
# A second elementary integration check of the same density series.
u=P.gen();cp=C*(2+u)**2;sp=S*(2+u)**2;fp=F*(4-u)**2
low=sum(r[j]/(2*j+4)+s[j]/(2*j+4)**2 for j in range(len(r)))
v=low+sum(2*cp[j]/(j+1) for j in range(0,N,2))+sum((sp[j]+fp[j])/(R(j)+R(3)/2) for j in range(N))
compare('four_squared_radius',mpnum(v),4,115)
for n in (3,4):
compare('normalization_%d'%n,vals[n][0],1,120)
matrix=mp.matrix([[vals[n][i+j] for j in range(5)] for i in range(5)])
mp.cholesky(matrix)
assert all(vals[n][k]>0 for k in (2,4,6,8))
evidence['checks']['moment_matrix_%d'%n]={'positive_definite':True,'orders':'0 through 8'}
# Numerically integrate the piecewise density, independently of the generator's
# coefficient products and elementary log-moment antiderivatives. This shares
# the density recurrence, so it is an arithmetic check, not a second density.
with mp.workdps(55):
R,P,N,r,s,C,S,F=density4_series(45)
@lru_cache(None)
def density_for_quad(x):
if x==0 or x==4:return mp.mpf(0)
return mpnum(series_density(x))
for k in range(9):
q=mp.quad(lambda x: density_for_quad(x)*mp.log(x)**k if x else mp.mpf(0),[0,1,2,3,4])
compare('four_numeric_quadrature_%d'%k,q,vals[4][k],40)
# Check readback against the independent numerical references retained here.
if '--stored' in sys.argv or os.environ.get('CHECK_STORED')=='1':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':numberdb.configure(api_key=sys.stdin.read().strip())
doc=numberdb.table('T686');rows=doc['Numbers'];count=0
assert set(rows)=={'3','4'}
stored_values={n:[mp.mpf(1),vals[n][1]] for n in (3,4)}
for n in (3,4):
for k in range(2,9):
assert set(rows[str(n)])==set(str(j) for j in range(2,9))
row=rows[str(n)][str(k)]
text=row['number'] if isinstance(row,dict) else row
stored_values[n].append(mp.mpf(text))
from decimal import Decimal
ulp=mp.power(10,Decimal(text).as_tuple().exponent)
error=abs(mp.mpf(text)-vals[n][k])
assert error<=ulp,(n,k,'readback exceeds one written unit')
evidence['checks']['stored_%d_%d'%(n,k)]={'error_in_written_ulps':mp.nstr(error/ulp,12)}
count+=1
assert count==14
for n in (3,4):
mp.cholesky(mp.matrix([[stored_values[n][i+j] for j in range(5)] for i in range(5)]))
evidence['checks']['stored_moment_matrix_%d'%n]={'positive_definite':True}
compare('stored_four_closed_form',stored_values[4][2],closed4,99)
for k in range(2,9):compare('stored_three_radial_%d'%k,stored_values[3][k],radial3(k),99)
report=HigherMahler().verify(sample=None);assert report.ok
evidence['verify']=str(report)
import sage.version
evidence['versions']={'mpmath':mp.__version__,'numberdb':numberdb.__version__,'SageMath':sage.version.version,'Python':sys.version}
evidence['zero_error_note']='0.0 means equality after rounding to the comparison precision, not exact equality of the mathematical numbers.'
print('EVIDENCE_JSON',json.dumps(evidence,sort_keys=True),flush=True)