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

Run with SageMath (or a suitable passagemath installation):
    sage -pip install numberdb
    sage -python generate.py              # verify the stored draft
    sage -python generate.py --fill       # fill with NUMBERDB_API_KEY set

The --fill operation writes entries to the existing draft. It does not offer,
review, or make the table public. Private reads also need an authorized key.

Janson (2007), equations (7), (194), (199), fixes the variance-t, unit-duration
normalization. All coefficients are exact rationals. Even moments are exact;
odd moments multiply a rational by an Arb enclosure of sqrt(pi/2).
The default selection is the first twelve positive moment orders; see
argument-selection.md. Working precision uses 180 guard bits; checks.py
measures the available precision and independently checks every stored value.
"""
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

WORKING_GUARD = 180
MAX_ORDER = 12

@lru_cache(maxsize=None)
def auxiliary(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, or M_n/sqrt(pi/2) for odd n, exactly."""
    if n < 0:
        raise ValueError('moment order must be nonnegative')
    q = auxiliary(n)[1][n]
    m = n//2
    if n%2:
        return QQ(factorial(n))*q/QQ(2**m*factorial(3*m+1))
    return QQ(factorial(n)*4**(3*m)*factorial(3*m))*q/QQ(2**m*factorial(6*m))

class MeanderMoments(numberdb.Generator):
    table = 'T526'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
             'derivation.md', 'skill-manifest.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 moment orders must be positive')
        factor = rational_factor(n)
        if n%2 == 0:
            return factor
        field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
        return field(factor)*(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()
    generator = MeanderMoments()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill twelve independently checked Brownian meander area moments'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)