back to table · edit · history · where entries came from · files · download
5114 bytes, as of the version from 2026-11-02 11:04 (current). Recorded here, not run.
"""Reciprocal-cubic Mahler measures -- numberdb.org/T658.
Run with SageMath and mpmath:
sage -pip install numberdb mpmath
sage -python generate.py # verify every stored entry
sage -python generate.py --publish # replace draft entries with API key
The publish method writes entries; it does not publish the draft to the public.
The sole parameter k is real. Its exact input identities are constructed anew
at each precision. The 37 selected arguments and their four mathematical roles
are explained in argument-selection.md. No result cache is needed to regenerate.
Jensen's formula reduces the normalized torus integral to 2/pi times the
integral of J_k(theta), 0 <= theta <= pi/2. With c=cos(theta), b=|k-4c^2|,
J=log(2c) if b<=4c, and J=log((b+sqrt(b^2-16c^2))/2) otherwise.
The latter expression has the finite limit log|k| at c=0, for k != 0.
All crossing points b=4c are inserted as endpoints. Tanh-sinh quadrature is
performed independently at 140 and 190 decimal digits, requiring its numerical
error estimates to be below 10^(-working+15) on every subinterval. These
estimates and precision agreement are not certified error bounds. The table
retains 100 significant digits with rigour heuristic (agreement-checked).
At k=0 the factorization (1+x)(1+y)(x+y) proves g(0)=0 exactly.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
import mpmath as mp
WORKING_DIGITS = (140, 190)
ALGEBRAIC = ('4+2*sqrt(5)', '4-2*sqrt(5)',
'8+4*sqrt(5)', '8-4*sqrt(5)')
def identities():
yield '0'
for n in range(1, 17):
yield str(n)
yield str(-n)
yield from ALGEBRAIC
def argument(identity):
if identity in ALGEBRAIC:
a, b = {'4+2*sqrt(5)': (4, 2), '4-2*sqrt(5)': (4, -2),
'8+4*sqrt(5)': (8, 4), '8-4*sqrt(5)': (8, -4)}[identity]
return mp.mpf(a) + b * mp.sqrt(5)
return mp.mpf(int(identity))
def quad_piece(f, a, b, working, depth=0):
value, error = mp.quad(f, [a, b], method='tanh-sinh',
error=True, maxdegree=9)
if mp.isfinite(value) and error < mp.power(10, -working + 15):
return value
if depth >= 6:
raise ArithmeticError('Quadrature failed to converge on a subinterval')
middle = (a + b) / 2
return (quad_piece(f, a, middle, working, depth + 1)
+ quad_piece(f, middle, b, working, depth + 1))
@lru_cache(maxsize=None)
def compute(identity, working):
with mp.workdps(working):
k = argument(identity)
if k == 0:
return '0'
edges = [mp.mpf(0), mp.pi / 2]
if -1 <= k < 0:
s = mp.sqrt(1 + k)
edges.extend(mp.acos(c) for c in ((1-s)/2, (1+s)/2))
elif 0 < k < 8:
edges.append(mp.acos((mp.sqrt(1+k)-1)/2))
edges = sorted(set(edges))
result = mp.mpf(0)
for left, right in zip(edges, edges[1:]):
c = mp.cos((left + right) / 2)
active = abs(k - 4*c*c) >= 4*c
def integrand(theta):
c = mp.cos(theta)
if not active:
return mp.log(2*c)
b = abs(k - 4*c*c)
# Rounding at an exact crossing can make D very slightly
# negative. This clips only endpoint roundoff on this branch.
discriminant = max(mp.mpf(0), b*b - 16*c*c)
return mp.log((b + mp.sqrt(discriminant)) / 2)
result += quad_piece(integrand, left, right, working)
return mp.nstr(2*result/mp.pi, working)
class ReciprocalCubic(numberdb.Generator):
table = 'T658'
parameters = ('k',)
type = 'R'
digits = 100
rigour = 'heuristic (agreement-checked)'
def enumerate(self):
for identity in identities():
yield {'k': identity}
def value(self, params, digits):
identity = str(params['k'])
if identity == '0':
return {'number': 0, 'comment': r'$P_0=(1+x)(1+y)(x+y)$, so $g(0)=0$ exactly.'}
result = {'number': numberdb.agreeing(
lambda working: compute(identity, working), at=WORKING_DIGITS)}
if identity in ALGEBRAIC:
result['param-latex'] = '$' + identity.replace('*sqrt(5)', r'\sqrt5') + '$'
if identity == '-1':
result['comment'] = r'$P_{-1}=(x+y+1)(xy+x+y)$, so $g(-1)=2\mu_3$ CITE{formula-special}.'
elif identity == '8':
result['comment'] = r'The triangular-lattice spanning-tree entropy, equal to twice the honeycomb-lattice spanning-tree entropy CITE{formula-special}.'
return result
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator = ReciprocalCubic()
if '--publish' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Compute ordinary reciprocal-cubic Mahler measures at 100 digits'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)