back to table · edit · history · where entries came from · files · download
6158 bytes, as of the version from 2026-11-02 11:01 (current). Recorded here, not run.
"""Independent checks for Hessian Mahler measures -- numberdb.org/T649.
Run: sage -python checks.py (with generate.py and values.json alongside).
Numerical checks use the original cubic, not the quadratic integrand used
by generate.py. At each x=e^(i theta), solve the cubic in y by Cardano's
formula and integrate sum(log(max(1,abs(y)))) over 0<=theta<=pi/3.
Symmetry gives the factor 3/pi. Split at all projected torus zeros and at
Cardano discriminant zeros; midpoint splits aid convergence. This check
is also floating quadrature, not certification.
"""
import json
import time
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
import mpmath as mp
from generate import ITEMS,mp_argument,exact_argument
def cubic_jensen(item,working=125):
mp.mp._tanh_sinh.clear()
with mp.workdps(working):
t=mp_argument(item)
k=mp.sign(t)*mp.root(abs(t),3)
omega=mp.exp(2*mp.pi*mp.j/3)
def integrand(theta):
x=mp.exp(mp.j*theta)
p=-k*x
q=x**3+1
if t==0:
return max(mp.mpf(0),mp.log(abs(q)))
if t==27:
ys=[-x-1,-omega*x-omega**2,-omega**2*x-omega]
return sum(max(mp.mpf(0),mp.log(abs(y))) for y in ys)
root=mp.sqrt(q*q/4+p**3/27)
choices=(-q/2+root,-q/2-root)
u=max(choices,key=abs)**(mp.mpf(1)/3)
v=-p/(3*u)
ys=[u+v,omega*u+omega**2*v,omega**2*u+omega*v]
return sum(max(mp.mpf(0),mp.log(abs(y))) for y in ys)
end=mp.pi/3
points=[mp.mpf(0),end]
if -1<t<27:
c=mp.acos((k-1)/2)
points.extend([c/3,min(2*c/3,(2*mp.pi-2*c)/3)])
if 0<t<27:
points.append(2*mp.acos(mp.sqrt(t/27))/3)
points=sorted(set(points))
# Remove coincident points within working-precision rounding.
unique=[points[0]]
for p in points[1:]:
if p-unique[-1]>mp.eps*100:
unique.append(p)
split=[unique[0]]
for a,b in zip(unique,unique[1:]):split.extend([(a+b)/2,b])
tolerance=mp.power(10,-working+15)
def adaptive(a,b,depth=0):
value,error=mp.quad(integrand,[a,b],error=True,maxdegree=8)
if error<tolerance*(b-a)/end:
return value,error
if depth>=24:raise ArithmeticError('Cubic check convergence: '+item['key'])
m=(a+b)/2
left,el=adaptive(a,m,depth+1);right,er=adaptive(m,b,depth+1)
return left+right,el+er
pieces=[adaptive(a,b) for a,b in zip(split,split[1:])]
result=sum(v for v,e in pieces);error=sum(e for v,e in pieces)
assert mp.isfinite(result)
return mp.nstr(result*3/mp.pi,working),mp.nstr(error*3/mp.pi,8)
def hypergeometric(item,working=125):
with mp.workdps(working):
t=mp_argument(item)
assert t<=-1 or t>=27
if abs(t)>27:
# Constant-term hypergeometric series. Its nth term is
# (3n)! / (3n*(n!)^3*t^n). The absolute ratio is <27/abs(t),
# so the remaining absolute tail is bounded geometrically.
# Explicit iteration avoids mpmath's expensive acceleration
# near the unit argument. Rounding remains non-certified.
term=2/t
total=term
n=1
ratio_bound=27/abs(t)
target=mp.power(10,-working+5)
while True:
term*=mp.mpf(3*n*(3*n+1)*(3*n+2))/((n+1)**3*t)
n+=1
total+=term
if abs(term)*ratio_bound/(1-ratio_bound)<target:break
if n>1000000:raise ArithmeticError('Series did not converge')
val=mp.log(abs(t))/3-total
else:
val=mp.log(abs(t))/3-mp.re(2/t*mp.hyper([mp.mpf(4)/3,mp.mpf(5)/3,1,1],[2,2,2],27/t))
return mp.nstr(val,working)
def agreement(a,b):
with mp.workdps(150):
delta=abs(mp.mpf(a)-mp.mpf(b))/max(1,abs(mp.mpf(b)))
return float(-mp.log10(delta)) if delta else 999.0
def symbolic_checks():
# c4^3/discriminant from the shared Weierstrass model.
R=PolynomialRing(QQ,'t');t=R.gen()
a2=-27*t**2;a4=216*t**3*(t-27);a6=-432*t**4*(t-27)**2
b2=4*a2;b4=2*a4;b6=4*a6;b8=4*a2*a6-a4*a4
c4=b2*b2-24*b4
disc=-b2*b2*b8-8*b4**3-27*b6**2+9*b2*b4*b6
assert c4**3*(t-27)**3==disc*t*(t+216)**3
# Torus monomial substitution, after clearing negative powers:
# (x^2*y)^3 P(y/x^2,1/(x*y)) = Q(x^3,y^3).
R3=PolynomialRing(QQ,('x','y','k'));x,y,k=R3.gens()
assert y**6+x**3+x**6*y**3-k*x**3*y**3 == y**6+(x**6-k*x**3)*y**3+x**3
def run(values):
symbolic_checks()
# A known, nonzero control for the independent integrator runs first.
with mp.workdps(130):
L=(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
control=mp.nstr(3*mp.sqrt(3)*L/(4*mp.pi),125)
zero=next(i for i in ITEMS if i['key']=='0')
c,err=cubic_jensen(zero)
assert agreement(c,control)>115
print(json.dumps({'control':'N(0) from Hurwitz zeta Dirichlet L-value','digits':agreement(c,control),'value':c}),flush=True)
for item in ITEMS:
key=item['key'];start=time.monotonic()
v,err=cubic_jensen(item)
digs=agreement(v,values[key]['high'])
assert digs>110,(key,digs,v,values[key]['high'])
row={'t':key,'cubic_jensen':v,'quadrature_estimate':err,'agreement_digits':digs}
if exact_argument(item)<=-1 or exact_argument(item)>27:
h=hypergeometric(item)
hd=agreement(h,values[key]['high'])
assert hd>110,(key,hd)
row.update(hypergeometric=h,hypergeometric_agreement_digits=hd)
if key=='27':
with mp.workdps(130):target=mp.nstr(3*mp.mpf(control),125)
assert agreement(v,target)>115
row['seconds']=round(time.monotonic()-start,2)
print(json.dumps(row),flush=True)
if __name__=='__main__':
values=json.loads(Path(__file__).with_name('values.json').read_text())
run(values)