"""Moments of the positive-part Brownian bridge area -- numberdb.org/T529.

Run with SageMath (or a suitable passagemath installation):
    sage -pip install numberdb
    sage -python generate.py             # verify all stored entries
    sage -python generate.py --fill      # fill an authorized draft

Set NUMBERDB_API_KEY to read or fill a private draft. In the campaign,
NUMBERDB_KEY_FROM_STDIN=1 reads that key from stdin and NUMBERDB_PUBLISH=1
selects --fill. Filling entries does not offer or publish the table.

The centered bridge has covariance min(s,t)-st on [0,1], and the area is
integral max(b(t),0) dt. Only raw moments of orders 1 through 12 are selected.
The rational recurrences are Janson (2007), equations (7), (240), (245).
See argument-selection.md for the selection and derivation.md for checks.
At 100 digits the 180 guard bits give 513 working bits; checks.py measures
a minimum relative accuracy of 509 bits and verifies the final decimal interval, not just its center.
"""
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(order):
    """Exact K_n and P_n=D_n^+ through the requested order."""
    k = [QQ(-1)/QQ(2)]
    p = [QQ(1)]
    for n in range(1, order+1):
        k.append(QQ(3*n-4)/QQ(4)*k[n-1]
                 + sum((k[j]*k[n-j] for j in range(1,n)), QQ(0)))
        p.append(sum((k[j]*p[n-j] for j in range(1,n+1)), QQ(0)))
    return tuple(k), tuple(p)


def rational_factor(n):
    """M_n, with sqrt(pi/2) removed when n is odd."""
    if n < 0:
        raise ValueError('The moment order must be nonnegative')
    p = coefficients(n)[1][n]
    m = n//2
    if n % 2:
        return QQ(factorial(n))*p / QQ(2**m*factorial(3*m+1))
    return QQ(factorial(n)*4**(3*m)*factorial(3*m))*p / QQ(2**m*factorial(6*m))


class PositiveBridgeMoments(numberdb.Generator):
    table = 'T529'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
             'derivation.md', 'numerical-evidence.json', 'skill-manifest.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 moment orders are positive')
        q = rational_factor(n)
        if n % 2 == 0:
            return q
        r = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
        return r(q)*(r.pi()/r(2)).sqrt()


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    generator = PositiveBridgeMoments()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        result = generator.publish(message='fill twelve independently checked positive-bridge moments')
        print(result)
        print('run:', result.run)
    else:
        result = generator.verify(sample=None)
        print(result)
        sys.exit(0 if result.ok else 1)