back to table · edit · history · where entries came from · files · download
6235 bytes, as of the version from 2026-11-09 09:32 (current). Recorded here, not run.
"""Higher Mahler measures of 1+x_1+...+x_(n-1) -- numberdb.org/T686.
Run with SageMath:
sage -pip install numberdb mpmath # once
sage -python generate.py # verify all stored entries
sage -python generate.py --fill # fill the existing private draft
The generator computes raw signed logarithmic moments, with normalized Haar
measure, natural logarithms, and no factorial divisor. The default selection
is n=3,4 and k=2,...,8. See argument-selection.md and method.md.
All arithmetic in truncated series uses Arb. Series tails are not certified:
150- and 220-decimal working computations are compared, and the declared
rigour is heuristic (agreement-checked), not proven.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from math import factorial
WORKING_DIGITS = (150, 220)
MAX_ORDER = 8
def falling(x, j):
z = QQ(1)
for i in range(j):
z *= x-i
return z
def ode_coefficients(center, direction):
P = PolynomialRing(QQ, 'u'); u = P.gen(); x = center+direction*u
return [list(x**4-64), list(direction*(7*x**5-32*x**3+64*x)),
list(6*x**6-60*x**4), list(direction*(x**7-20*x**5+64*x**3))]
def frobenius(a, exponent, seeds, count, R):
"""Coefficient recurrence from BSWZ (11), expanded about 2 or 4."""
c = [R(v) for v in seeds]
r = QQ(exponent)
for m in range(len(c), count):
denom = R(0); rhs = R(0)
for j, poly in enumerate(a):
for ell, v in enumerate(poly):
if not v:
continue
i = m+j-ell-2
if i < 0 or i > m:
continue
z = v*falling(r+i, j)
if i == m:
denom += z
else:
rhs += z*c[i]
assert not denom.contains_zero()
c.append(-rhs/denom)
return c
def low4(count, R):
"""BSWZ (26),(30); p4(x)=sum (r_j-s_j log(x))*x^(2j+1)."""
s = [3/(2*R.pi()**2)]
r = [s[0]*R(8).log()]
for j in range(1, count):
a=4*(2*j-1)*(5*j*j-5*j+2); b=-2*(j-1)**3
sm2=s[j-2] if j>1 else R(0)
rm2=r[j-2] if j>1 else R(0)
sj=(a*s[-1]+b*sm2)/(128*j**3)
rj=(a*r[-1]+b*rm2+3*(64*j*j*sj-(20*j*j-20*j+6)*s[-1]+(j-1)**2*sm2))/(128*j**3)
s.append(sj);r.append(rj)
return r,s
def density4_series(working):
R=RealBallField(numberdb.bits(3*working+50, losing=96))
P=PolynomialRing(R,'u'); N=4*working
a=ode_coefficients(2,1)
A=P(frobenius(a,0,[1,0],N,R))
B=P(frobenius(a,0,[0,1],N,R))
S=P(frobenius(a,QQ(1)/2,[1],N,R))*(-4/R.pi()**2)
p2=R(2)**(QQ(7)/3)*R.pi()/(3*R(3).sqrt()*R(QQ(2)/3).gamma()**6)
p1=R(5).sqrt()/(40*R.pi()**4)
for j in (1,2,4,8):p1*=R(QQ(j)/15).gamma()
slope=(p1-p2*A(-1))/B(-1)
C=p2*A+slope*B
F=P(frobenius(ode_coefficients(4,-1),QQ(1)/2,[R(2).sqrt()/R.pi()**2],N,R))
r,s=low4(2*working,R)
return R,P,N,r,s,C,S,F
def moments4(working, max_order=MAX_ORDER):
R,P,N,r,s,C,S,F=density4_series(working)
log2=P([R(2).log()]+[R((-1)**(j+1))/(j*2**j) for j in range(1,N)])
log4=P([R(4).log()]+[-R(1)/(j*4**j) for j in range(1,N)])
pow2=P(1);pow4=P(1);values=[]
for k in range(max_order+1):
low=sum((r[j]*((-1)**k*factorial(k))/R(2*j+2)**(k+1)
-s[j]*((-1)**(k+1)*factorial(k+1))/R(2*j+2)**(k+2)) for j in range(len(r)))
cp=(C*pow2).truncate(N);sp=(S*pow2).truncate(N);fp=(F*pow4).truncate(N)
mid=sum(2*cp[j]/(j+1) for j in range(0,N,2))
high=sum((sp[j]+fp[j])/(R(j)+QQ(3)/2) for j in range(N))
values.append(low+mid+high)
pow2=(pow2*log2).truncate(N);pow4=(pow4*log4).truncate(N)
return values
def moments3(working,max_order=MAX_ORDER):
"""Taylor coefficients of BBSW Theorem 4.1, equation (23)."""
R=RealBallField(numberdb.bits(working,losing=96));P=PolynomialRing(R,'s');s=P.gen()
N=max_order+1;term=P(1);total=term
for j in range(2*working):
# (1+s/2+j)^3 / (4*(j+1)^2*(j+3/2+s/2)).
b=R(j)+QQ(3)/2
inverse=P([(-R(1)/(2*b))**k/b for k in range(N)])
term=(term*(j+1+s/2)**3*inverse).truncate(N)/(4*(j+1)**2)
total+=term
# log(3^(s+1) Gamma(1+s/2)^2/Gamma(2+s)).
logarithm=[R(3).log(),R(3).log()-1]
for k in range(2,N):
logarithm.append(((-1)**k*(R(2)**(1-k)-1)*R(k).zeta()+(-1)**k)/k)
exponential=[R(3)]
for k in range(1,N):
exponential.append(sum(j*logarithm[j]*exponential[k-j] for j in range(1,k+1))/k)
result=(P(exponential)*total).truncate(N)*R(3).sqrt()/(2*R.pi())
return [result[k]*factorial(k) for k in range(N)]
@lru_cache(None)
def computed(n, working):
v=(moments3 if n==3 else moments4)(working)
# These are decimal approximations to truncated computations. Their Arb
# radii bound arithmetic only, not the omitted series tails.
import mpmath as mp
with mp.workdps(working+10):
return tuple(mp.nstr(mp.mpf(x.mid().str(truncate=False)),working) for x in v)
class HigherMahler(numberdb.Generator):
table='T686'
parameters=('n','k')
type='R'
digits=100
rigour='heuristic (agreement-checked)'
def enumerate(self):
for n in (3,4):
for k in range(2,9):
yield {'n':n,'k':k}
def value(self, params, digits):
n,k=int(params['n']),int(params['k'])
if n not in (3,4) or not 2<=k<=8:
raise ValueError('Implemented selection: n=3,4 and 2<=k<=8.')
return numberdb.agreeing(lambda w:computed(n,w)[k],at=WORKING_DIGITS)
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator=HigherMahler()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
print(generator.publish(message='Computed 14 higher Mahler measures with independent controls'))
else:
report=generator.verify(sample=None);print(report)
sys.exit(0 if report.ok else 1)