"""Broadhurst–Mellit determinants of Bessel moments -- numberdb.org/T647.

Run with SageMath and the numberdb Python package:
    sage -pip install numberdb
    sage -python generate.py                 # verify all stored rows
    sage -python generate.py --fill          # write entries to this draft
The second command needs NUMBERDB_API_KEY. It fills entries; it does not
publish or offer a draft. NUMBERDB_KEY_FROM_STDIN=1 reads a key from stdin;
NUMBERDB_PUBLISH=1 selects the same fill operation for agents/sage.sh.

Zhou, arXiv:1711.01829v2, equations (1.4), (1.6), proved in Section 4.
For n >= 3 the determinant has order floor((n-1)/2), rows a=1,...,m,
and columns r=1,3,...,2m-1 of integral I_0(t)^a K_0(t)^(n-a) t^r dt.
All scales are 1. The determinant is signed, without added normalizations.

The only free parameter n is the discrete total factor count. The default
range 3..22 covers both parities at orders 1..10. Its finite endpoint is an
editorial reference choice, not a mathematical boundary; see the attached
argument-selection.md. There are no selection data files or imported local
helpers. Dependencies: Python standard library, numberdb, SageMath (Arb).

Every division of exact prefactors is over QQ. The even formula's Gamma
argument and exponent are rational, not Python floats. RealBallField keeps
rigorous enclosures throughout. At 100 requested digits, 64 extra bits give
397-bit working precision. checks.py measures retained accuracy and checks
all rows at a higher precision, plus independently encoded recursions.
"""
import os
import sys
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

WORKING_GUARD = 64
DEFAULT_MAX_N = 22


def determinant(n, digits=100, guard=WORKING_GUARD):
    n = ZZ(n)
    if n < 3:
        raise ValueError('n must be an integer at least 3')
    k = (n - 1) // 2
    R = RealBallField(numberdb.bits(digits, losing=guard))
    pi = R.pi()
    if n % 2:
        numerator = ZZ(1)
        radicand = ZZ(1)
        for j in range(1, int(k) + 1):
            numerator *= ZZ(2*j)**(k-j)
            radicand *= ZZ(2*j+1)**(2*j+1)
        result = R(numerator)*pi**(k*(k+1)//2)/R(radicand).sqrt()
    else:
        prefactor = QQ(2)
        for j in range(1, int(k) + 2):
            prefactor *= QQ(ZZ(2*j-1)**(k+1-j))/QQ(ZZ(2*j)**j)
        exponent = QQ((k+1)**2)/QQ(2)
        gamma_argument = QQ(k+1)/QQ(2)
        result = R(prefactor)*pi**exponent/R(gamma_argument).gamma()
    if not result.is_finite() or not result > 0:
        raise ArithmeticError('determinant must be a finite positive ball')
    return result


class BroadhurstMellit(numberdb.Generator):
    table = 'T647'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'proven'

    def enumerate(self, max_n=DEFAULT_MAX_N):
        for n in range(3, max_n + 1):
            yield {'n': n}

    def value(self, params, digits):
        return determinant(params['n'], digits)


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=', 1)[0].isidentifier():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    g = BroadhurstMellit()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(g.publish(message='Fill 20 certified full Bessel moment determinants'))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)