program.py

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

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

import numberdb.sage
from math import factorial
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

def moment(n):
    """Raw moment of the standard Brownian excursion area, n >= 1."""
    if n < 1:
        raise ValueError('n must be positive')
    k = [-QQ(1)/QQ(2)]
    for j in range(1, n+1):
        k.append(QQ(3*j-4)/QQ(4)*k[j-1]
                 + sum((k[i]*k[j-i] for i in range(1,j)), QQ(0)))
    m = n//2
    if n%2:
        r = RealBallField(512)
        q = QQ(4*factorial(n))*k[n]/QQ(2**m*factorial(3*m))
        return r(q)*(r.pi()/r(2)).sqrt()
    j = 3*m-1
    return QQ(4*factorial(n)*4**j*factorial(j))*k[n]/QQ(2**m*factorial(2*j))

# For example, the mean:
print(moment(1))