generate.py

back to table · edit · history · where entries came from · files · download

3892 bytes, as of the version from 2026-11-02 00:20. Recorded here, not run.

"""Mahler measures of x+x^-1+y+y^-1+z+z^-1-k -- numberdb.org/T663.

Run with SageMath:
    sage -pip install numberdb mpmath
    sage -python generate.py                 # verify every stored row
    sage -python generate.py --fill          # fill this existing draft

--fill requires NUMBERDB_API_KEY; it writes values, never changes draft status.
The two working precisions are 150 and 210 decimal digits. Quadrature agreement
is heuristic, not a certified enclosure. See argument-selection.md for the
14 exact integer arguments and numerical-evidence.md for measured checks.

Jensen integrates out z. The sum s=2cos(u)+2cos(v) has probability density
K(1-s^2/16)/(2*pi^2) = 1/(4*pi*AGM(1,abs(s)/4)), for -4<s<4,
where K uses the elliptic parameter. Only abs(k-s)>2 contributes. Splitting
at s=0 and at the endpoints k+-2 handles logarithmic and square-root behavior.
The AGM spelling avoids rounding 1-s^2/16 to 1 near the logarithmic singularity.
"""
import sys
import numberdb.sage as numberdb
import mpmath as mp

ARGUMENTS = tuple(range(13)) + (18,)
WORKING_DIGITS = (150, 210)


def density(s):
    # A value assigned at an isolated integrable singularity changes no integral.
    if s == 0:
        return mp.mpf(0)
    return 1 / (4 * mp.pi * mp.agm(1, abs(s) / 4))


def compute(k, working):
    """Ordinary logarithmic measure for a finite real coefficient k."""
    with mp.workdps(working):
        k = abs(mp.mpf(k))
        if not mp.isfinite(k):
            raise ValueError('A finite real coefficient is required.')
        total = mp.mpf(0)
        pieces = [(-mp.mpf(4), min(mp.mpf(4), k - 2), 1),
                  (max(-mp.mpf(4), k + 2), mp.mpf(4), -1)]
        for left, right, sign in pieces:
            if left >= right:
                continue
            edges = [left] + ([mp.mpf(0)] if left < 0 < right else []) + [right]
            def integrand(s):
                t = sign * (k - s) - 2
                # acosh(1+t/2)=2asinh(sqrt(t)/2), stable at t=0.
                return density(s) * 2 * mp.asinh(mp.sqrt(max(0, t)) / 2)
            value, error = mp.quad(integrand, edges, error=True, maxdegree=10)
            if not mp.isfinite(value) or error > mp.power(10, -working + 20):
                raise ArithmeticError('Quadrature convergence estimate failed.')
            total += value
        return mp.nstr(total, working)


COMMENTS = {
    0: "Smyth's value $m(P_0)=m(1+x+y)$ CITE{formula-zero}.",
    2: r'The interior critical value of the torus sum; it is also an integer CM example in CITE{BF}, Section 2.3.',
    3: r'An integer CM example CITE{BF}, Section 2.3.',
    6: r'The simple-cubic spanning-tree entropy CITE{formula-lattice}; this is the boundary of the torus-zero interval.',
    10: r'An integer CM example CITE{BF}, Section 2.3.',
    18: r'An integer CM example CITE{BF}, Section 2.3.',
}


class BertinP(numberdb.Generator):
    table = 'T663'
    files = ('generate.py', 'checks.py', 'argument-selection.md',
             'numerical-evidence.md', 'methods.md')
    parameters = ('k',)
    type = 'R'
    digits = 100
    rigour = 'heuristic (agreement-checked)'

    def enumerate(self):
        for k in ARGUMENTS:
            yield {'k': str(k)}

    def value(self, params, digits):
        if digits > 100:
            raise ValueError('This measured protocol supports 100 stored digits.')
        k = int(params['k'])
        value = numberdb.agreeing(lambda working: compute(k, working), at=WORKING_DIGITS)
        result = {'number': value}
        if k in COMMENTS:
            result['comment'] = COMMENTS[k]
        return result


if __name__ == '__main__':
    g = BertinP()
    if '--fill' in sys.argv:
        print(g.publish(message='Fill ordinary Bertin P Mahler measures after independent checks'))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)