"""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)