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

WORKING_GUARD = 180

@lru_cache(maxsize=None)
def coefficients(n):
    k, q = [-QQ(1)/QQ(2)], [QQ(1)]
    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)))
        q.append(QQ(3*j-2)/QQ(2)*q[j-1]
                 + 2*sum((k[i]*q[j-i] for i in range(1,j+1)), QQ(0)))
    return tuple(k), tuple(q)

def rational_factor(n):
    """Return M_n for even n, or M_n / sqrt(2/pi) for odd n."""
    _, q = coefficients(n)
    w = sum((q[j]*q[n-j] for j in range(n+1)), QQ(0))
    m = n//2
    if n % 2 == 0:
        return QQ(factorial(n))*w / QQ(2**m*factorial(3*m))
    j = 3*m+2
    return QQ(factorial(n)*4**j*factorial(j))*w / QQ(2**(m+1)*factorial(2*j))

def moment(n, digits=100):
    if not isinstance(n, int) or n < 0:
        raise ValueError('moment order must be a nonnegative Python integer')
    a = rational_factor(n)
    if n % 2 == 0:
        return a
    r = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    return r(a)*(r(2)/r.pi()).sqrt()

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

The recurrences for $K_n$ and $Q_n$, the convolution for $W_n$, and the rational prefactors were evaluated exactly over $\mathbb Q$. Even moments are stored as exact rational numbers. Odd moments use an Arb enclosure of $\sqrt{2/\pi}$ at 513 bits, with a measured minimum relative accuracy of 509 bits, and are written to 100 significant digits. Each decimal denotes the interval within one unit of its last place.