"""Independent gamma/elliptic checks for T521's five even additions.

    python even-checks.py even-computation.jsonl

No generator or Bessel evaluator is imported. BBBG, arXiv:0801.0891v2,
equations (47), (49), (78), (89), (95), and (100) give the identities used.
The quarter-interval reduction immediately following (81) removes the
logarithmic endpoint from c40. Combining (89) and (95) gives
c42=(-2*c40+3*B1+3*pi^2-3*H)/16. Elliptic K takes the parameter m here,
whereas the paper writes its modulus. Carlson RF avoids subtracting tiny
complementary parameters in the s40 integral. Independent results are
agreement-checked at 140 and 180 digits, not labelled certified enclosures.
"""

import json
import sys
from decimal import Decimal
from pathlib import Path

import mpmath as mp


EVEN_KEYS = ((0, 3, 0), (0, 3, 2), (0, 4, 0), (0, 4, 2), (1, 3, 0))


def reference_values(working):
    with mp.workdps(working):
        third_gamma = mp.gamma(mp.mpf(1) / 3)
        third_zero = 3 * third_gamma ** 6 / (32 * mp.pi * mp.root(4, 3))
        third_second = third_zero / 9 - mp.pi ** 4 / (24 * third_zero)

        def elliptic_square(angle):
            return mp.ellipk(mp.sin(angle) ** 2) ** 2

        fourth_zero = 4 * mp.pi * mp.quad(
            elliptic_square, [0, mp.pi / 8, mp.pi / 4], method='gauss-legendre')

        def subtracted(angle):
            if not angle:
                return mp.pi ** 3 / 2
            return (4 * mp.pi * elliptic_square(angle) - mp.pi ** 3) / mp.sin(angle) ** 2

        subtracted_integral = mp.quad(
            subtracted, [0, mp.pi / 8, mp.pi / 4], method='gauss-legendre')
        gamma_constant = mp.pi ** 3 - (
            mp.gamma(mp.mpf(1) / 4) ** 4 + 16 * mp.gamma(mp.mpf(3) / 4) ** 4) / 8
        fourth_second = (-2 * fourth_zero + 3 * subtracted_integral
                         + 3 * mp.pi ** 2 - 3 * gamma_constant) / 16

        def complementary_product(variable):
            if not variable:
                return mp.mpf(0)
            angle = variable ** 2
            sine_squared = mp.sin(angle) ** 2
            cosine_squared = mp.cos(angle) ** 2
            return (4 * variable * mp.elliprf(0, cosine_squared, 1)
                    * mp.elliprf(0, sine_squared, 1))

        mixed_zero = mp.quad(complementary_product,
                            [0, mp.sqrt(mp.pi / 8), mp.sqrt(mp.pi / 4)],
                            method='tanh-sinh')
        result = dict(zip(EVEN_KEYS, (third_zero, third_second, fourth_zero,
                                     fourth_second, mixed_zero)))
        if not all(mp.isfinite(value) and value > 0 for value in result.values()):
            raise ArithmeticError('Nonfinite or nonpositive independent value')
        return result


def check(rows):
    entries = {tuple(int(row['params'][name]) for name in ('a', 'b', 'k')): row
               for row in rows}
    if not set(EVEN_KEYS).issubset(entries):
        raise ValueError('Missing even moments')
    lower = reference_values(140)
    higher = reference_values(180)
    records = []
    with mp.workdps(210):
        for key in EVEN_KEYS:
            reference = higher[key]
            agreement = abs(lower[key] - reference) / abs(reference)
            produced = mp.mpf(entries[key]['at_200'])
            difference = abs(produced - reference) / abs(reference)
            written = entries[key]['number']
            text_center = mp.mpf(written)
            decimal_exponent = Decimal(written).as_tuple().exponent
            written_error_ulps = abs(text_center - reference) / mp.mpf(10) ** decimal_exponent
            if not agreement < mp.mpf('1e-125') or not difference < mp.mpf('1e-125'):
                raise ArithmeticError(('Independent check failed', key, agreement, difference))
            if not written_error_ulps < mp.mpf('0.6'):
                raise ArithmeticError(('Stored last-place claim failed', key, written_error_ulps))
            records.append({
                'params': key, 'independent_at_140': mp.nstr(lower[key], 140),
                'independent_at_180': mp.nstr(reference, 180),
                'independent_precision_agreement': mp.nstr(agreement, 12),
                'relative_difference_from_bessel_quadrature': mp.nstr(difference, 12),
                'stored_error_ulps': mp.nstr(written_error_ulps, 12),
                'finite_and_positive': True})
    return records


if __name__ == '__main__':
    rows = [json.loads(line) for line in Path(sys.argv[1]).read_text().splitlines() if line]
    print(json.dumps({'independent_checks': check(rows), 'complete': True}, indent=2))