back to table · edit · history · where entries came from · files · download
8893 bytes, as of the version from 2026-11-02 11:02 (current). Recorded here, not run.
"""Ordinary and modified Hessian periods at exact real algebraic arguments.
Run with SageMath and the numberdb package installed. Each argument is given
by an integer polynomial and an exact rational isolating interval. It is
reconstructed at each working precision, never from a stored decimal.
Jensen's formula reduces the ordinary measure to the integral of the
logarithm of the larger modulus of the roots of
y^2 + (x^2-k*x)*y + x, where k is the real cube root of t.
The modified period is I-2J, with the split at acos((k-1)/2).
Precision agreement and adaptive quadrature estimates are not certified
error bounds. The value at t=0 of the modified period is exactly zero.
"""
from functools import lru_cache
import math
import numberdb.sage as numberdb
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ
from sage.rings.qqbar import AA
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.real_mpfr import RealField
import mpmath as mp
LOW_EXTRA_DIGITS = 40
HIGH_EXTRA_DIGITS = 90
MINIMUM_AGREEMENT_MARGIN = 15
INTEGRATION_EXTRA_DIGITS = 40
@lru_cache(maxsize=512)
def algebraic_argument(coefficients, interval):
polynomial_ring = PolynomialRing(QQ, 'parameter')
polynomial = polynomial_ring(list(coefficients))
lower, upper = map(QQ, interval)
if lower == upper:
if polynomial(lower) != 0:
raise ValueError('The specified rational argument is not a root')
return AA(lower)
if lower >= upper:
raise ValueError('Invalid isolating interval')
bounds = RealIntervalField(256)(lower, upper)
result = AA.polynomial_root(polynomial, bounds)
if not lower < result < upper:
raise ValueError('The selected root is outside its isolating interval')
return result
def exact_argument(item):
return algebraic_argument(tuple(item['polynomial']), tuple(item['isolating_interval']))
def numerical_argument(item, working):
field = RealField(math.ceil((working + 40) * math.log2(10)) + 64)
return mp.mpf(str(field(exact_argument(item))))
def three_step_measure():
third = mp.mpf(1) / 3
dirichlet_value = (mp.zeta(2, third) - mp.zeta(2, 2 * third)) / 9
return 3 * mp.sqrt(3) * dirichlet_value / (4 * mp.pi)
def jensen_integrand(coefficient, angle):
sine = mp.sin(angle / 2)
difference = coefficient - 3
imaginary = sine * (difference + 4 * sine ** 2) / 2
modulus_difference = (difference * (coefficient + 1) + 4 * coefficient * sine ** 2) / 4
if imaginary == 0:
return mp.asinh(mp.sqrt(max(modulus_difference, mp.mpf(0))))
discriminant = mp.sqrt(modulus_difference ** 2 + 4 * imaginary ** 2)
if modulus_difference < 0:
hyperbolic_sine = abs(imaginary) * mp.sqrt(2 / (discriminant - modulus_difference))
else:
hyperbolic_sine = mp.sqrt((discriminant + modulus_difference) / 2)
return mp.asinh(hyperbolic_sine)
def compute_periods(item, working, force_integral=False):
mp.mp._tanh_sinh.clear()
exact = exact_argument(item)
with mp.workdps(working + INTEGRATION_EXTRA_DIGITS):
argument = numerical_argument(item, working + INTEGRATION_EXTRA_DIGITS)
interior = bool(-1 < exact < 27)
if exact == 0 and not force_integral:
return {'ordinary': mp.nstr(three_step_measure(), working - 5),
'modified': '0', 'method': 'Dirichlet control and exact cancellation',
'evaluations': 0, 'quadrature_estimate': None}
if exact == 27 and not force_integral:
return {'ordinary': mp.nstr(3 * three_step_measure(), working - 5),
'modified': None, 'method': 'singular-fibre factorization',
'evaluations': 0, 'quadrature_estimate': None}
coefficient = mp.sign(argument) * mp.root(abs(argument), 3)
evaluations = 0
transition_points = []
if coefficient:
sine_squared = (3 - coefficient) * (coefficient + 1) / (4 * coefficient)
if 0 < sine_squared < 1:
center = 2 * mp.asin(mp.sqrt(sine_squared))
sine = mp.sin(center / 2)
imaginary = sine * (coefficient - 3 + 4 * sine ** 2) / 2
scale = 4 * abs(imaginary) / (abs(coefficient) * mp.sin(center))
if scale:
transition_points.append(center)
while scale < mp.pi:
transition_points.extend(point for point in (center - scale, center + scale)
if 0 < point < mp.pi)
scale *= 4
def integrand(angle):
nonlocal evaluations
evaluations += 1
result = jensen_integrand(coefficient, angle)
if not mp.isfinite(result):
raise ArithmeticError('Nonfinite integrand')
return result
tolerance = mp.power(10, -working + 20)
def integrate(lower, upper, depth=0):
if lower == upper:
return mp.mpf(0), mp.mpf(0)
value, estimate = mp.quad(integrand, [lower, upper], error=True, maxdegree=8)
if mp.isfinite(value) and estimate <= tolerance * (upper - lower) / mp.pi:
return value, estimate
if depth >= 24:
raise ArithmeticError('Quadrature exhausted its subdivision limit for ' + item['key'])
middle = (lower + upper) / 2
first, first_estimate = integrate(lower, middle, depth + 1)
second, second_estimate = integrate(middle, upper, depth + 1)
return first + second, first_estimate + second_estimate
def interval(lower, upper):
points = sorted({lower, upper, (lower + upper) / 2}
| {point for point in transition_points if lower < point < upper})
values = [integrate(left, right) for left, right in zip(points, points[1:])]
return mp.fsum(value for value, estimate in values) / mp.pi, mp.fsum(estimate for value, estimate in values) / mp.pi
if interior:
crossing = mp.acos((coefficient - 1) / 2)
first, first_estimate = interval(mp.mpf(0), crossing)
second, second_estimate = interval(crossing, mp.pi)
ordinary = first + second
modified = first - 2 * second
estimate = first_estimate + second_estimate
else:
ordinary, estimate = interval(mp.mpf(0), mp.pi)
modified = None
if not mp.isfinite(ordinary) or ordinary <= 0:
raise ArithmeticError('Invalid ordinary Mahler measure')
if modified is not None and not mp.isfinite(modified):
raise ArithmeticError('Invalid modified Mahler measure')
if estimate > tolerance:
raise ArithmeticError('Insufficient quadrature convergence')
result = {
'ordinary': mp.nstr(ordinary, working - 5),
'modified': mp.nstr(modified, working - 5) if modified is not None else None,
'method': 'quadratic Jensen integral', 'evaluations': evaluations,
'quadrature_estimate': mp.nstr(estimate, 12),
'integration_extra_digits': INTEGRATION_EXTRA_DIGITS,
'transition_breakpoints': len(transition_points),
}
mp.mp._tanh_sinh.clear()
return result
def evaluate(item, digits=100):
precisions = (digits + LOW_EXTRA_DIGITS, digits + HIGH_EXTRA_DIGITS)
measurements = [compute_periods(item, precision) for precision in precisions]
result = {'key': item['key'], 'digits': digits, 'working_precisions': list(precisions),
'measurements': measurements, 'values': {}, 'agreement_digits': {}}
with mp.workdps(digits + HIGH_EXTRA_DIGITS + 30):
for quantity in ('ordinary', 'modified'):
if measurements[1][quantity] is None:
continue
if quantity == 'modified' and exact_argument(item) == 0:
result['values'][quantity] = 0
result['agreement_digits'][quantity] = 'exact'
continue
lower, higher = [mp.mpf(measurement[quantity]) for measurement in measurements]
if not mp.isfinite(lower) or not mp.isfinite(higher) or higher == 0:
raise ArithmeticError('Nonfinite or unsupported zero result')
discrepancy = abs((lower - higher) / higher)
agreement = -mp.log10(discrepancy) if discrepancy else mp.inf
if agreement < digits + MINIMUM_AGREEMENT_MARGIN:
raise ArithmeticError('Insufficient precision agreement for ' + item['key'] + ': ' + str(agreement))
result['agreement_digits'][quantity] = str(agreement)
result['values'][quantity] = mp.nstr(higher, digits, strip_zeros=False)
return result