"""Verify the right-triangle index clarification for numberdb.org/T485.

Run with SageMath after installing numberdb:
    sage -python check_t485_repair.py T485-document.json
For a private live draft, pipe its authorized API key on standard input:
    cat "$NUMBERDB_KEY_FILE" | sage -python check_t485_repair.py

Read-only: independently enumerate the two sorted integer-index multisets,
including multiplicities, and enclose their first 20 values at 512 bits.
No attached generator is imported. Decimal intervals use NumberDB's one-unit
last-place convention; zeros must be exact. The j cutoff is checked against
the last included coefficient, so omitted indices cannot enter the first 20.
References: Band--Bersudsky--Fajman (arXiv:1507.03410v1), Section 2;
Laugesen--Pan--Son (doi:10.7153/mia-15-32), Section 6 and Table 2.
"""
import json
import sys
import urllib.request
from decimal import Decimal
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION

if len(sys.argv) > 1:
    with open(sys.argv[1]) as source:
        document = json.load(source)
else:
    key = sys.stdin.read().strip()
    request = urllib.request.Request(
        'https://numberdb.org/api/table?id=T485',
        headers={'Authorization': 'Bearer ' + key})
    with urllib.request.urlopen(request) as response:
        document = json.load(response)

R = RealBallField(512)
pi_squared = R.pi() ** 2
bound = 30
checks = []
for domain in ('right-isosceles-triangle', 'hemiequilateral-triangle'):
    hemi = domain == 'hemiequilateral-triangle'
    scale = QQ(16) / QQ(9) if hemi else QQ(1)
    modes = sorted((scale * (j*j + k*k + (j*k if hemi else 0)), j, k)
                   for j in range(bound + 1) for k in range(j + 1))
    stored = document['Numbers'][domain]
    assert list(stored) == [str(n) for n in range(1, 21)]
    assert modes[0][0] == 0
    assert modes[1][0] == scale  # first nonzero coefficient control
    assert scale * (bound + 1)**2 > modes[19][0]
    for n, (coefficient, j, k) in enumerate(modes[:20], 1):
        value = stored[str(n)]['number']
        if coefficient == 0:
            assert value == '0'
        else:
            dec = Decimal(value)
            midpoint = Fraction(dec)
            unit = Fraction(10) ** dec.as_tuple().exponent
            lo, hi = midpoint - unit, midpoint + unit
            lower, upper = QQ(lo.numerator)/QQ(lo.denominator), QQ(hi.numerator)/QQ(hi.denominator)
            ball = R(coefficient) * pi_squared
            assert ball.is_finite()
            assert ball > R(lower) and ball < R(upper), (domain, n)
        checks.append({'domain': domain, 'n': n, 'j': j, 'k': k,
                       'coefficient_of_pi_squared': str(coefficient),
                       'stored_interval_contains_ball': True})

print(json.dumps({'sage_version': SAGE_VERSION, 'working_bits': 512,
                  'checked': len(checks), 'exact_zeros': 2,
                  'tail_bound_verified': True, 'rows': checks}, indent=2))