"""Mehta integrals for Gaussian ensembles -- numberdb.org/T566.

Run with SageMath (or a compatible passagemath installation):
    sage -pip install numberdb
    sage -python generate.py                    # verify every stored row
    sage -python generate.py --fill             # write entries to this draft
    sage -python checks.py                      # independent integral checks

--fill requires NUMBERDB_API_KEY. It only writes entries to the existing
private draft; it does not offer, review or publish the table. To pass a key
through a wrapper, set NUMBERDB_KEY_FROM_STDIN=1 and pipe the key on stdin;
NUMBERDB_PUBLISH=1 selects --fill for wrappers that do not forward arguments.

The full labelled-coordinate integral has weight exp(-sum(x_i^2)/2).
Gamma factors are reduced exactly to q * sqrt(2)^s * sqrt(pi)^p before any
ball arithmetic. For n>=1 the power p is positive, so no selected row is
rational. The selection is beta=1,2,4 and n=1,...,50; see argument-selection.md.
A 96-bit guard is measured by checks.py, which also verifies the serialized
100-digit intervals against independent exact Gaussian moment integrals.
"""
import os
import sys
from math import factorial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

MAX_N = 50
DYSON_EXPONENTS = (1, 2, 4)
WORKING_GUARD = 96


def half_gamma(twice_argument):
    """Return (q,p) for Gamma(k/2)=q*sqrt(pi)^p, k a positive integer."""
    k = int(twice_argument)
    if k < 1:
        raise ValueError('Gamma argument must be positive')
    if k % 2 == 0:
        return QQ(factorial(k // 2 - 1)), 0
    m = (k - 1) // 2
    return QQ(factorial(2 * m)) / QQ(4 ** m * factorial(m)), 1


def exact_components(beta, n):
    """Return q,s,p for Z=q*sqrt(2)^s*sqrt(pi)^p with s=0 or 1."""
    beta, n = int(beta), int(n)
    if beta not in DYSON_EXPONENTS or n < 1:
        raise ValueError('Expected beta in (1,2,4) and n>=1')
    q = QQ(2) ** (n // 2)
    p = n
    denominator, dp = half_gamma(2 + beta)
    for j in range(1, n + 1):
        numerator, np = half_gamma(2 + j * beta)
        q *= numerator / denominator
        p += np - dp
    assert q > 0 and p > 0
    return q, n % 2, p


def evaluate_components(parts, precision):
    q, s, p = parts
    field = RealBallField(precision)
    result = field(q) * field(2).sqrt() ** s * field.pi().sqrt() ** p
    if not result.is_finite() or not result > 0:
        raise ArithmeticError('Nonpositive or nonfinite integral enclosure')
    return result


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

    def enumerate(self, max_n=MAX_N):
        for beta in DYSON_EXPONENTS:
            for n in range(1, int(max_n) + 1):
                yield {'beta': str(beta), 'n': str(n)}

    def value(self, params, digits):
        parts = exact_components(params['beta'], params['n'])
        result = evaluate_components(parts, numberdb.bits(digits, losing=WORKING_GUARD))
        if int(result.accuracy()) < numberdb.bits(digits):
            raise ArithmeticError('Insufficient relative accuracy')
        return result


def authenticated_client():
    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('\"\'')
        return numberdb.Client(api_key=key)
    return numberdb.Client()


if __name__ == '__main__':
    generator = MehtaIntegrals()
    client = authenticated_client()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        outcome = generator.publish(client=client, message='Fill 150 Mehta integrals certified by independent Gaussian moment checks')
        print(outcome)
    report = generator.verify(sample=None, client=client)
    print(report)
    sys.exit(0 if report.ok else 1)