generate.py

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

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

"""Moments of the absolute Brownian bridge area -- numberdb.org/T524.

Run with SageMath (or a modular Sage installation):
    sage -pip install numberdb
    sage -python generate.py             # verify every stored row
    sage -python generate.py --fill      # fill this draft using NUMBERDB_API_KEY

For a private draft, verification also needs an authorized API key.
--fill writes entries, not publication or review state.
The only argument is positive integer raw moment order; the default is 1..12.
See argument-selection.md for the standard unit-duration bridge normalization.
Janson (2007), doi:10.1214/07-PS104, equations (135) and (149), supplies
an exact rational recurrence. Half-integer Gamma factors cancel before
conversion: even moments are exact rationals and odd moments are rational
multiples of sqrt(pi/2). Arb supplies certified enclosures for the latter.
The 180-bit guard gives 513 working bits at 100 digits; checks.py measures
relative accuracy and validates the actual serialized intervals.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from functools import lru_cache
from math import factorial
import os
import sys

MAX_ORDER = 12
WORKING_GUARD = 180


@lru_cache(maxsize=None)
def d_coefficients(maximum):
    d = [QQ(1)]
    for n in range(1, maximum + 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))
    return tuple(d)


def rational_factor(n):
    """M_n for even n, or M_n/sqrt(pi/2) for odd n; also accepts seed 0."""
    if n < 0:
        raise ValueError('moment order must be nonnegative')
    d = d_coefficients(n)[n]
    m = n // 2
    if n % 2:
        return QQ(factorial(n))*d/QQ(2**m*factorial(3*m+1))
    return QQ(factorial(n)*4**(3*m)*factorial(3*m))*d/QQ(2**m*factorial(6*m))


class AbsoluteBridgeAreaMoments(numberdb.Generator):
    table = 'T524'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
             'numerical-evidence.json', 'skill-read-evidence.json')

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

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


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    g = AbsoluteBridgeAreaMoments()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(g.publish(message='Fill twelve standard Brownian bridge absolute-area moments'))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)