generate.py

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

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

"""Moments of the Brownian excursion area -- numberdb.org/T523.

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

The raw moments use standard Brownian variance t and duration 1. The only
parameter is positive integer moment order; enumerate() selects 1 through 12.
See argument-selection.md for the selection and its source.

Janson (2007), doi:10.1214/07-PS104, equations (6)-(7), supplies the rational
K recurrence. Half-integer Gamma factors are cancelled exactly. Even moments
are returned in QQ; odd moments are rational multiples of sqrt(pi/2), computed
in Arb with at least 512 bits. At 100 digits the guard is 180 bits beyond
numberdb.bits(100). checks.py measures actual relative accuracy and validates
the 100-significant-digit written enclosures. Filling does not publish or offer
the draft. The executable check and selection note are attached with this file.
"""
import os
import sys
from functools import lru_cache
from math import factorial

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

MAX_ORDER = 12
WORKING_GUARD = 180


@lru_cache(maxsize=None)
def takacs(n):
    """Exact K_n from Janson (6)-(7)."""
    if n == 0:
        return -QQ(1) / QQ(2)
    return (QQ(3*n - 4) / QQ(4) * takacs(n - 1)
            + sum((takacs(j) * takacs(n - j) for j in range(1, n)), QQ(0)))


def rational_factor(n):
    """Return M_n for even n, or M_n/sqrt(pi/2) for odd n."""
    if n < 1:
        if n == 0:
            return QQ(1)
        raise ValueError('moment order must be nonnegative')
    m = n // 2
    if n % 2:
        return QQ(4 * factorial(n)) * takacs(n) / QQ(2**m * factorial(3*m))
    j = 3*m - 1
    return (QQ(4 * factorial(n) * 4**j * factorial(j)) * takacs(n)
            / QQ(2**m * factorial(2*j)))


class ExcursionMoments(numberdb.Generator):
    table = 'T523'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md')

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

    def value(self, params, digits):
        n = int(params['n'])
        if n < 1:
            raise ValueError('stored moment order must be positive')
        q = rational_factor(n)
        if n % 2 == 0:
            return q
        field = RealBallField(max(512, numberdb.bits(digits, losing=WORKING_GUARD)))
        result = field(q) * (field.pi() / field(2)).sqrt()
        assert result.is_finite() and result > 0
        return result


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=', 1)[0].strip().isidentifier():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    generator = ExcursionMoments()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Twelve raw moments: exact rationals and certified 100-digit enclosures'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)