"""Bessel moments, T521: the original quadrature extended to five even moments.

    python -m pip install numberdb python-flint mpmath
    python generate.py --compute --new-only
    python generate.py --compute
    python generate.py

The default enumeration contains all 40 original odd moments and five even
moments specified in even-selection.md. The last command verifies the live
table without writing. No command publishes, reviews, or approves anything.

This extends the attached generator with SHA256
6f50e2a959efd804b670232fe50fdb6ab25048e388278cc83652438a94aa9cfc:
the same direct tanh-sinh quadrature, Arb Bessel evaluations at exact binary
nodes, 150/200 working digits, radius gates, and exponential tail estimate.
Python-flint replaces Sage's Arb interface; no special-function method is
replaced. Old decimal strings are retained only after reproduction checks.

For k=0 only, t=u^2 changes the integrand to
2*u*I0(u^2)^a*K0(u^2)^b, whose limit at zero is zero. This avoids falsely
assigning a finite limiting value to K0(t)^b itself. Positive k uses the
original variable and panels. Both cases truncate at the same t=T.

For t>=T>=1, I0(t)<=exp(t), K0(t)<=sqrt(pi/(2*t))*exp(-t), b<=6,
b-a>=1 and 0<=k<=7 imply the common bound
(pi/2)^3 * 7! * exp(-T) * sum(T^j/j!, j=0..7)
for the omitted tail. These inequalities follow from DLMF 10.32.1 and
10.32.9, using cosh(v)>=1+v^2/2 for K0. T=3*(working+60).
The substitution does not change the integral or its tail at t=T.

Arb midpoints enter mpmath quadrature; discretization is NOT certified.
The honest label remains heuristic (agreement-checked), at 100 digits.
Independent gamma and elliptic-integral identities are in even-checks.py.
"""

import argparse
import json
from functools import lru_cache

import mpmath as mp
import numberdb
from flint import arb, ctx


WORKING_DIGITS = (150, 200)
STORED_DIGITS = 100
EVEN_KEYS = ((0, 3, 0), (0, 3, 2), (0, 4, 0), (0, 4, 2), (1, 3, 0))


def selection():
    for total in range(3, 7):
        for first in range((total - 1) // 2 + 1):
            second = total - first
            powers = {1, 3, 5, 7}
            powers.update(power for left, right, power in EVEN_KEYS
                          if (left, right) == (first, second))
            for power in sorted(powers):
                yield {'a': first, 'b': second, 'k': power}


def binary_midpoint(value):
    mantissa, exponent = value.man_exp()
    return mp.ldexp(mp.mpf(int(mantissa)), int(exponent))


@lru_cache(maxsize=65536)
def _bessels(binary_node, working):
    sign, mantissa, exponent, _ = binary_node
    node = mp.mp.make_mpf(binary_node)
    previous_precision = ctx.prec
    ctx.prec = numberdb.bits(working + 40) + 3 * int(mp.ceil(node))
    try:
        exact_node = arb((-1 if sign else 1) * mantissa) * arb(2) ** exponent
        values = []
        for ball in (exact_node.bessel_i(0), exact_node.bessel_k(0)):
            if not ball.is_finite():
                raise ArithmeticError('Nonfinite Bessel enclosure')
            midpoint = binary_midpoint(ball.mid())
            radius = binary_midpoint(ball.rad())
            if not mp.isfinite(midpoint) or midpoint <= 0:
                raise ArithmeticError('Bessel midpoint is not finite and positive')
            if not radius < midpoint * mp.mpf(10) ** (-working - 15):
                raise ArithmeticError('Insufficient Bessel enclosure accuracy')
            values.append(midpoint)
        return tuple(values)
    finally:
        ctx.prec = previous_precision


def tail_bound(working):
    with mp.workdps(working):
        endpoint = 3 * (working + 60)
        return ((mp.pi / 2) ** 3 * mp.factorial(7) * mp.exp(-endpoint)
                * mp.fsum(mp.mpf(endpoint) ** power / mp.factorial(power)
                          for power in range(8)))


def integral(params, working):
    first, second, power = (int(params[name]) for name in ('a', 'b', 'k'))
    if not (0 <= first < second <= 6 and 0 <= power <= 7):
        raise ValueError('This bound covers 0<=a<b<=6 and 0<=k<=7 only')
    with mp.workdps(working):
        def integrand(variable):
            if not variable:
                return mp.mpf(0)
            node = variable ** 2 if power == 0 else variable
            jacobian = 2 * variable if power == 0 else 1
            first_value, second_value = _bessels(node._mpf_, working)
            return jacobian * node ** power * first_value ** first * second_value ** second

        endpoint = 3 * (working + 60)
        panels = [0, 1, 4, 16, 64, endpoint]
        if power == 0:
            panels = [mp.sqrt(point) for point in panels]
        value, estimated_error = mp.quad(integrand, panels, method='tanh-sinh', error=True)
        if not mp.isfinite(value) or value <= 0:
            raise ArithmeticError('Nonfinite or nonpositive moment')
        if estimated_error >= abs(value) * mp.mpf(10) ** (-working + 10):
            raise ArithmeticError('Quadrature convergence estimate failed')
        if tail_bound(working) >= abs(value) * mp.mpf(10) ** (-working):
            raise ArithmeticError('Tail bound too large')
        return +value


def evaluated(params, digits=STORED_DIGITS):
    if digits > STORED_DIGITS:
        raise ValueError('This run justifies at most 100 stored digits')
    low = integral(params, WORKING_DIGITS[0])
    high = integral(params, WORKING_DIGITS[1])
    with mp.workdps(WORKING_DIGITS[1]):
        relative = abs(low - high) / abs(high)
        if relative >= mp.mpf('1e-130'):
            raise ArithmeticError('150/200-digit agreement failed')
        return {
            'params': params, 'number': mp.nstr(high, digits, strip_zeros=False),
            'at_150': mp.nstr(low, 150), 'at_200': mp.nstr(high, 200),
            'relative_difference': mp.nstr(relative, 12),
            'upper_tail_bounds': [mp.nstr(tail_bound(working), 12)
                                  for working in WORKING_DIGITS],
            'endpoint_substitution': 't=u^2' if int(params['k']) == 0 else 'none',
            'rigour': 'heuristic (agreement-checked)'}


class BesselMoments(numberdb.Generator):
    table = 'T521'
    parameters = ('a', 'b', 'k')
    type = 'R'
    digits = STORED_DIGITS
    rigour = 'heuristic (agreement-checked)'
    files = ('generate.py', 'even-checks.py', 'even-selection.md')

    def enumerate(self):
        yield from selection()

    def value(self, params, digits):
        return evaluated(params, digits)['number']


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--compute', action='store_true')
    parser.add_argument('--new-only', action='store_true')
    options = parser.parse_args()
    if options.compute:
        for parameters in selection():
            identity = tuple(parameters[name] for name in ('a', 'b', 'k'))
            if not options.new_only or identity in EVEN_KEYS:
                print(json.dumps(evaluated(parameters)), flush=True)
    else:
        if options.new_only:
            parser.error('--new-only requires --compute')
        report = BesselMoments().verify(sample=None)
        print(report)
        raise SystemExit(0 if report.ok else 1)