"""Moments of the Brownian double meander area -- numberdb.org/T527.

Run with SageMath (the numberdb package must be installed):
    sage -pip install numberdb
    sage -python generate.py              # verify all stored entries
    sage -python generate.py --fill       # fill this draft with NUMBERDB_API_KEY

The family uses standard Brownian motion of duration one, minus its GLOBAL
minimum. Only positive integer raw moment orders 1..12 are enumerated.
Janson (2007), equations (7), (199), (229), (232), gives the recurrences.
Exact QQ arithmetic simplifies Gamma at integers and half-integers. Odd
moments are rational multiples of sqrt(2/pi), enclosed by Arb. The precision
is bits(100, losing=180)=513 bits; measured accuracy is recorded in the
attached numerical-evidence.json. Every decimal is checked as its implied
interval, not just at its midpoint. See argument-selection.md and checks.py.
"""
import os
import sys
from math import factorial
from functools import lru_cache
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 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()

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

    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 orders must be positive')
        value = moment(n, digits)
        if n == 1:
            return {'number': value,
                    'comment': r'The mean is $\sqrt{2/\pi}$, the expected maximum of standard Brownian motion on $[0,1]$.'}
        return value

if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if key.startswith('NUMBERDB_API_KEY='):
            key = key.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    gen = DoubleMeanderMoments()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(gen.publish(message='Compute twelve raw moments with exact rational factors and Arb enclosures'))
    else:
        report = gen.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)