generate.py

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

3494 bytes, as of the version from 2026-10-26 19:49 (current). Recorded here, not run.

"""Moments of the absolute Brownian motion area -- numberdb.org/T525.

Run with SageMath (or passagemath with numberdb installed):
    sage -pip install numberdb
    sage -python generate.py             # verify every stored entry
    sage -python checks.py               # offline independent checks
    sage -python generate.py --fill      # write entries to the existing draft

--fill uses NUMBERDB_API_KEY; it does not publish or offer the table.
For a private draft, verification also requires its author's API key.
NUMBERDB_KEY_FROM_STDIN=1 reads a token from stdin for campaign wrapper use.

Janson, Probability Surveys 4 (2007), equations (149), (181), (190),
arXiv:0704.2289, DOI 10.1214/07-PS104. D and L are computed over QQ.
Integer/half-integer Gamma factors are simplified exactly. Even moments
are rational; odd moments are q*sqrt(2/pi) enclosed with Arb. Working
precision is bits(digits, losing=180), 513 bits at 100 requested digits.
The guard is measured by checks.py, including the final decimal intervals.
The generator has no selection-dependent imports or input files.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from math import factorial
import os
import sys

WORKING_GUARD = 180
MAX_ORDER = 12


def coefficients(limit):
    """Return exact D and L from the two cited rational recurrences."""
    d, ell = [QQ(1)], [QQ(1)]
    for n in range(1, limit + 1):
        d.append(QQ(3*n-2)/QQ(4)*d[n-1]
                 - sum((d[j]*d[n-j] for j in range(1, n)), QQ(0))/QQ(2))
        ell.append(QQ(3*n-1)/QQ(2)*ell[n-1]
                   - sum((d[n-j]*ell[j] for j in range(1, n)), QQ(0)))
    return d, ell


def rational_factor(n, ell):
    """M_n itself for even n; M_n/sqrt(2/pi) for odd n."""
    m = n // 2
    if n % 2 == 0:
        return QQ(factorial(n))*ell / QQ(2**m * factorial(3*m))
    j = 3*m + 2
    return QQ(factorial(n)*4**j*factorial(j))*ell / QQ(2**(m+1)*factorial(2*j))


def moment(n, digits=100):
    n = int(n)
    if n < 0:
        raise ValueError('moment order must be nonnegative')
    ell = coefficients(n)[1][n]
    q = rational_factor(n, ell)
    if n % 2 == 0:
        return q
    r = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    return r(q) * (r(2)/r.pi()).sqrt()


class AbsoluteBrownianAreaMoments(numberdb.Generator):
    table = 'T525'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
             'derivation.md', 'skill-manifest.json')

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

    def value(self, params, digits):
        n = int(params['n'])
        if n < 1:
            raise ValueError('stored moment orders must be positive')
        return moment(n, digits)


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 = AbsoluteBrownianAreaMoments()
    if '--fill' in sys.argv:
        print(generator.publish(message='Fill twelve independently checked Brownian absolute-area moments'))
    report = generator.verify(sample=None)
    print(report)
    sys.exit(0 if report.ok else 1)