generate.py

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

3719 bytes, as of the version from 2026-11-09 16:36 (current). Recorded here, not run.

"""Taylor coefficients of the Bettin-Conrey reciprocity function -- numberdb.org/T692.

Run with SageMath and the numberdb Python package:
    sage -pip install numberdb             # once
    sage -python generate.py               # verify all stored entries
    sage -python generate.py --publish     # submit entries to this draft
The second command requires NUMBERDB_API_KEY. It does not offer or publish a
website draft. Reading an unpublished draft also requires its author's key.
NUMBERDB_KEY_FROM_STDIN=1 reads a private key from stdin; the
campaign wrapper also accepts NUMBERDB_PUBLISH=1 for the second command.

AM (1.3), https://arxiv.org/html/2402.14216v1, gives the finite expression.
Euler's even-zeta formula makes each coefficient an exact polynomial in pi^2.
This implementation uses exact Sage Bernoulli numbers and rational polynomial
coefficients, followed by Arb Horner evaluation. All divisions are exact Sage
rational operations. 384 extra bits cover the measured cancellation through
order 100. Each returned ball must retain at least the requested digits plus
32 extra bits; this assertion also guards extensions. The two rational values
are known from the definition and reciprocal symmetry, not inferred from digits.
"""
import os
import sys
from functools import lru_cache
from math import comb, factorial
import numberdb.sage as numberdb
from sage.arith.misc import bernoulli
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

MAX_ORDER = 100
WORKING_GUARD = 384


@lru_cache(None)
def rational_coefficients(n):
    """Coefficients of P_n(X) with g_n=P_n(pi^2), in ascending degree."""
    if n < 0:
        raise ValueError('Taylor order must be nonnegative')
    if n == 0:
        return (QQ(-1),)
    if n == 1:
        return (QQ(1)/2,)
    p = [QQ(0) for _ in range(n//2+1)]
    p[0] = QQ(1)/(n*(n+1))
    for k in range(2,n+1,2):
        # b_k=B_k*zeta(k)/k=(-1)^(k/2+1)*2^(k-1)*B_k^2*pi^k/(k*k!).
        b = QQ((-1)**(k//2+1)*2**(k-1))*bernoulli(k)**2/(k*factorial(k))
        weight = 2*comb(n-1,k-2) + (2 if k == n else 0)
        p[k//2] += weight*b
    return tuple(p)


def coefficient(n, digits=100):
    p = rational_coefficients(n)
    if n <= 1:
        return p[0]
    R = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    x = R.pi()**2
    result = R(0)
    for c in reversed(p):
        result = result*x+R(c)
    assert result.is_finite()
    assert result.accuracy() >= numberdb.bits(digits, losing=32)
    return result


class BettinConreyTaylor(numberdb.Generator):
    table = 'T692'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'proven'

    def enumerate(self):
        for n in range(MAX_ORDER+1):
            yield {'n': n}

    def value(self, params, digits):
        n = int(params['n'])
        v = coefficient(n,digits)
        if n == 0:
            return {'number': v, 'comment': r'$g_0=-1$, since $c_0(1)=0$.'}
        if n == 1:
            return {'number': v, 'comment': r'$g_1=1/2$, by $g(x)=xg(1/x)$ at $x=1$.'}
        return v


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if key.startswith(('NUMBERDB_API_KEY=','NUMBERDB_KEY=')):
            key = key.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    generator = BettinConreyTaylor()
    if '--publish' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill the normalized Taylor coefficient prefix with certified values'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)