back to table · edit · history · where entries came from · files · download
6002 bytes, as of the version from 2026-11-02 16:53. Recorded here, not run.
"""Higher Mahler measure of cyclotomic polynomials -- numberdb.org/T681.
Run with SageMath (or a compatible passagemath installation):
sage -pip install numberdb
sage -python generate.py # verify every stored entry
sage -python generate.py --fill # fill an existing draft with a key
For a private draft set NUMBERDB_API_KEY, or set NUMBERDB_KEY_FROM_STDIN=1
and pipe the key into the command. --fill sends entries; it does not publish
or offer the draft. The default selection is every odd squarefree n <= 105,
then k=2,3. See the attached argument-selection.md for the coverage rationale.
Lalin--Sinha, arXiv:1106.1304v1, Theorems 5 and 6 supply the finite formulas.
Their degree is called d here; n always means the cyclotomic order. The exact
primitive-root angles are a/n, not rounded angles from a table. C3 and S2
mean Re Li_3(exp(i theta)) and Im Li_2(exp(i theta)). Reflection is applied
before evaluating the functions. Rational pair weights are combined exactly.
WORKING_GUARD=128 supplies extra precision for the cancellation in Theorem 6.
The full default range retained at least 414 relative accuracy bits
(about 124 decimal digits) at 461 working bits in this build. Every result must retain the requested precision plus 32 bits; if a
larger future range needs more, the computation increases precision explicitly.
All transcendental operations are Arb enclosures, not wrapped approximations.
"""
import os
import sys
from collections import defaultdict
from functools import lru_cache
from math import gcd
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
WORKING_GUARD = 128
MAX_ORDER = 105
MOMENT_ORDERS = (2, 3)
def factor_integer(n):
"""Positive integer factorization; no floating arithmetic."""
result = []
p = 2
while p * p <= n:
exponent = 0
while n % p == 0:
exponent += 1
n //= p
if exponent:
result.append((p, exponent))
p += 1
if n > 1:
result.append((n, 1))
return result
def orders(limit=MAX_ORDER):
return [n for n in range(1, int(limit) + 1, 2)
if all(e == 1 for _, e in factor_integer(n))]
def second_coefficient(n):
"""The exact coefficient of pi^2 in LS Theorem 5."""
coefficient = QQ(1) / 12
for p, _ in factor_integer(n):
coefficient *= QQ(2 * (p - 1)) / p
return coefficient
@lru_cache(maxsize=None)
def pair_weights(n):
residues = [a for a in range(n) if gcd(a, n) == 1]
degree = len(residues)
cosine = defaultdict(int)
sine = defaultdict(lambda: QQ(0))
for i, a in enumerate(residues):
for j in range(i + 1, degree):
delta = QQ(residues[j] - a) / n
weight = degree * delta - (j - i)
if delta > QQ(1) / 2:
delta = 1 - delta
weight = -weight
cosine[delta] += 1
sine[delta] += weight
return degree, dict(cosine), dict(sine)
@lru_cache(maxsize=16384)
def clausen_pair(delta, precision):
R = RealBallField(precision)
C = ComplexBallField(precision)
if delta == 0:
return R(3).zeta(), R(0)
if delta == QQ(1) / 2:
return -R(3).zeta() * QQ(3) / 4, R(0)
z = (C(0, 1) * 2 * C.pi() * C(delta)).exp()
return z.polylog(3).real(), z.polylog(2).imag()
def moment_at_precision(n, k, precision):
R = RealBallField(precision)
if k == 2:
return R.pi() ** 2 * second_coefficient(n)
if k != 3:
raise ValueError('This generator implements only moment orders 2 and 3')
degree, cosine_weights, sine_weights = pair_weights(n)
csum = R(0)
ssum = R(0)
for delta, multiplicity in cosine_weights.items():
c3, s2 = clausen_pair(delta, precision)
csum += multiplicity * c3
ssum += sine_weights[delta] * s2
return (-QQ(3) / 2 * degree ** 2 * R(3).zeta()
- 3 * degree * csum - 3 * R.pi() * ssum)
def moment(n, k, digits=100):
n, k = int(ZZ(n)), int(ZZ(k))
if n < 1:
raise ValueError('Cyclotomic order n must be positive')
precision = numberdb.bits(digits, losing=WORKING_GUARD)
required = numberdb.bits(digits, losing=32)
for _ in range(5):
result = moment_at_precision(n, k, precision)
if result.is_finite() and result.accuracy() >= required:
return result
precision += 64
raise ArithmeticError('The finite formula did not reach the requested accuracy')
class HigherCyclotomicMahler(numberdb.Generator):
table = 'T681'
parameters = ('n', 'k')
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'quadrature.py', 'check.py', 'flatten.py',
'argument-selection.md')
def enumerate(self, limit=MAX_ORDER):
for n in orders(limit):
for k in MOMENT_ORDERS:
yield {'n': n, 'k': k}
def value(self, params, digits):
n, k = int(params['n']), int(params['k'])
value = moment(n, k, digits)
if n == 1:
return {'number': value, 'comment': (
r'$\Phi_1(x)=x-1$; the moment is $\pi^2/12$.' if k == 2 else
r'$\Phi_1(x)=x-1$; the moment is $-3\zeta(3)/2$.')}
return value
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
token = sys.stdin.read().strip()
if '=' in token and token.split('=', 1)[0].isidentifier():
token = token.split('=', 1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY'] = token
generator = HigherCyclotomicMahler()
if '--fill' in sys.argv:
print(generator.publish(message='Compute signed higher Mahler moments in Arb'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)