back to table · edit · history · where entries came from · files · download
3301 bytes, as of the version from 2026-11-02 17:11 (current). Recorded here, not run.
"""Independent direct-integral checks for T681; no writes.
Run: sage -python quadrature.py, with generate.py alongside it.
Mpmath tanh-sinh quadrature uses the Mobius product, not either finite
Clausen formula. This is a heuristic numerical check, NOT the certification
of the stored digits. The generator's finite Arb formulas certify those.
The integral is split at every j/n. All logarithms of sine are evaluated
with exact small residue numerators and endpoint distances, avoiding a loss
from evaluating sin(pi*m) close to an integer m. Removable zeros are split
as well, so no quadrature interval crosses a root of a Mobius factor.
"""
import json
from math import gcd
import time
import mpmath as mp
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.arith.misc import moebius
from numberdb._write import to_text
from generate import moment
def integral(n, dps):
with mp.workdps(dps):
terms=[(n//int(d),int(moebius(n//d))) for d in ZZ(n).divisors()
if moebius(n//d)]
total=mp.mpc(0)
for j in range((n+1)//2):
def integrand(u):
if not u or u==1:
return mp.mpc(0)
f=mp.mpf(0)
for q,mu in terms:
y=(j%q)+u
if y>mp.mpf(q)/2:
# Preserve 1-u, rather than cancelling at a large q.
y=(q-1-j%q)+(1-u)
f+=mu*mp.log(2*mp.sinpi(y/q))
return mp.mpc(f*f,f*f*f)
top=mp.mpf(1)/2 if j==(n-1)//2 else mp.mpf(1)
total += mp.quad(integrand,[0,top],method='tanh-sinh')
return total*2/n
def main():
ns=[n for n in range(1,106,2) if moebius(n)!=0]
# Nonzero positive and negative known controls must work first.
mp.mp.dps=60
generating = lambda s: mp.gamma(1+s)/mp.gamma(1+s/2)**2
gamma_second=mp.diff(generating,0,2)
gamma_third=mp.diff(generating,0,3)
assert abs(gamma_second-mp.pi**2/12)<mp.mpf('1e-55')
assert abs(gamma_third+mp.mpf(3)*mp.zeta(3)/2)<mp.mpf('1e-55')
ctrl=integral(1,50)
assert abs(ctrl.real-mp.pi**2/12)<mp.mpf('1e-42')
assert abs(ctrl.imag+mp.mpf(3)*mp.zeta(3)/2)<mp.mpf('1e-42')
print('CONTROL: direct n=1 integral matches pi^2/12 and -3*zeta(3)/2 to 42 decimal places',flush=True)
rows=[]
for n in ns:
started=time.monotonic()
val=ctrl if n==1 else integral(n,50)
targets=[mp.mpf(to_text(moment(n,k),digits=100)) for k in (2,3)]
errors=[abs(val.real-targets[0]),abs(val.imag-targets[1])]
assert max(errors)<mp.mpf('1e-40'), (n,errors)
row={'n':n,'m2':mp.nstr(val.real,45),'m3':mp.nstr(val.imag,45),
'absolute_errors':[mp.nstr(e,8) for e in errors],
'seconds':round(time.monotonic()-started,3)}
rows.append(row); print('QUAD',json.dumps(row),flush=True)
print('QUADRATURE_JSON '+json.dumps({'table':'T681','dps':50,'absolute_tolerance':'1e-40',
'method':'mpmath tanh-sinh; Mobius log-sine sum; split at j/n; combined second and third moments',
'mpmath_version':mp.__version__,'values':rows,
'skill_manifest':'agents/runs/20261102T162707Z-build-skills/manifest.json'}),flush=True)
if __name__=='__main__': main()