back to table · edit · history · where entries came from · files · download
3553 bytes, as of the version from 2026-10-26 22:20 (current). Recorded here, not run.
"""Moments of the positive-part Brownian motion area -- numberdb.org/T530.
Run with SageMath or passagemath and the numberdb client:
sage -pip install numberdb
sage -python generate.py # verify every stored entry
sage -python checks.py # independent local checks
sage -python generate.py --fill # fill this draft, API key required
The model and enumeration are justified in argument-selection.md. Janson
(2007), equations (7), (199), (245), (298), and (266), supply the exact rational
recurrences and moment formula. Rational half-integer Gamma factors are
simplified before evaluation. Only sqrt(2/pi) needs Arb. Working precision
is numberdb.bits(digits, losing=180), giving 513 bits for 100 digits; the
measured minimum retained relative accuracy is recorded in numerical-evidence.json.
The --fill option writes entries to this existing draft; it never publishes
or offers the table. NUMBERDB_KEY_FROM_STDIN=1 reads the key from stdin.
"""
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
WORKING_GUARD = 180
MAX_ORDER = 12
@lru_cache(maxsize=None)
def coefficients(last):
"""Return K, Q, P=D^+, and L^+ over QQ through order last."""
k = [QQ(-1)/2]
q = [QQ(1)]
p = [QQ(1)]
for n in range(1, last + 1):
k.append(QQ(3*n-4)/4*k[n-1]
+ sum(k[j]*k[n-j] for j in range(1, n)))
q.append(QQ(3*n-2)/2*q[n-1]
+ 2*sum(k[j]*q[n-j] for j in range(1, n+1)))
p.append(sum(k[j]*p[n-j] for j in range(1, n+1)))
ell = [(p[n] + sum(p[j]*q[n-j] for j in range(n+1)))/QQ(2)
for n in range(last+1)]
return tuple(k), tuple(q), tuple(p), tuple(ell)
def rational_factor(n):
"""M_n itself for even n, and M_n / sqrt(2/pi) for odd n."""
if n < 0:
raise ValueError('Moment order must be nonnegative')
ell = coefficients(n)[3][n]
if n % 2 == 0:
return QQ(factorial(n))*ell / (QQ(2)**(n//2)*factorial(3*n//2))
m = (n-1)//2
h = 3*m+2
return (QQ(factorial(n))*ell*QQ(4)**h*factorial(h)
/ (QQ(2)**(m+1)*factorial(2*h)))
class PositiveBrownianMoments(numberdb.Generator):
table = 'T530'
parameters = ('n',)
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'check.py', 'one_value.py', 'argument-selection.md',
'derivation.md', 'numerical-evidence.json', '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 start at one')
factor = rational_factor(n)
if n % 2 == 0:
return factor
field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
return field(factor)*(field(2)/field.pi()).sqrt()
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
generator = PositiveBrownianMoments()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill twelve raw moments with exact rational and Arb arithmetic'))
else:
result = generator.verify(sample=None)
print(result)
sys.exit(0 if result.ok else 1)