"""Check the amplitude normalization -- numberdb.org/T482.

Run with SageMath (no writes):
    sage -pip install numberdb
    sage -python normalization_checks.py
For a private draft, set NUMBERDB_API_KEY, or pipe the key on stdin with
NUMBERDB_KEY_FROM_STDIN=1. Output is a JSON check report.

The coefficient identity
 [x^n](1 - sqrt(1-4*x))/2 = Catalan(n-1)
is checked exactly for n=1,...,40. Stirling's formula then gives
 Catalan(n-1) ~ 4^n n^(-3/2)/(4*sqrt(pi)).
Thus the amplitude includes the square-root transfer factor 1/(2*sqrt(pi));
it is not the singular coefficient b, which is 1/2 in this example.
The general factor is -1/Gamma(-1/2) = 1/(2*sqrt(pi)).

The four closed-form entries are independently computed at 640 bits and
compared to the full stored decimal intervals (one unit in the last place).
Finite, narrow balls and an exact coefficient control prevent vacuous tests.
This does not re-certify the ten implicit-function amplitude enclosures.
"""
import json
import math
import os
import sys
import urllib.request
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.version import version as sage_version

B = RealBallField(640)
key = (sys.stdin.read().strip() if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1'
       else os.environ.get('NUMBERDB_API_KEY', ''))
headers = {'Authorization': 'Bearer ' + key} if key else {}
with urllib.request.urlopen(urllib.request.Request(
        'https://numberdb.org/api/table?id=T482', headers=headers), timeout=60) as r:
    doc = json.load(r)
assert 'Numbers' in doc, 'T482 is unavailable; private drafts require an API key.'

def rational(s):
    f = Fraction(s)
    return QQ(f.numerator) / QQ(f.denominator)

# Independent exact coefficient check for the generating function.
coefficient = QQ(1)
for n in range(1, 41):
    coefficient *= (QQ(1)/2 - (n-1)) * QQ(-4) / QQ(n)
    catalan = QQ(math.comb(2*(n-1), n-1)) / QQ(n)
    assert -coefficient/2 == catalan

transfer = -1 / B(-QQ(1)/2).gamma()
expected_transfer = 1 / (2 * B.pi().sqrt())
assert transfer.is_finite() and expected_transfer.is_finite()
assert B(transfer.rad()) < B(10)**-185 and B(expected_transfer.rad()) < B(10)**-185
assert (transfer - expected_transfer).contains_zero()

forms = {
    'plane': 1 / (4 * B.pi().sqrt()),
    'plane-binary': B(3).sqrt() / (2 * B.pi().sqrt()),
    'labelled-rooted': 1 / (2 * B.pi()).sqrt(),
    'labelled-free': 1 / (2 * B.pi()).sqrt(),
}
checks = {}
for name, value in forms.items():
    row = next(iter(doc['Numbers'][name]['nodes'].values()))
    s = row['number']
    assert '.' in s and 'e' not in s.lower()
    center = rational(s)
    ulp = QQ(1) / QQ(10)**len(s.split('.')[1])
    assert value.is_finite() and B(value.rad()) < B(ulp)*B(10)**-80
    assert B(center-ulp) < value and value < B(center+ulp)
    # A missing transfer factor must fail the same stored interval test.
    if name == 'plane':
        assert B(center+ulp) < B(QQ(1)/2)
    checks[name] = {'contained_in_stored_last_place_interval': True,
                   'midpoint_discrepancy_ulps': str(abs(value-B(center))/B(ulp)),
                   'computed_ball': str(value)}
ratios = []
for n in (200, 400, 800):
    count = QQ(math.comb(2*(n-1), n-1)) / QQ(n)
    normalized = B(count) * B(n)**(QQ(3)/2) / B(4)**n
    ratio = normalized/forms['plane']
    assert ratio.is_finite() and B(1) < ratio < B(QQ(1003)/1000)
    ratios.append({'nodes': n, 'normalized_count_over_amplitude': str(ratio)})
assert len(doc['Numbers']) == 14
assert all(len(v['nodes']) == 1 for v in doc['Numbers'].values())
print(json.dumps({'sage_version': sage_version, 'precision_bits': 640,
                  'exact_Catalan_coefficient_checks': 40,
                  'transfer_factor_identity': True,
                  'closed_form_entries': checks,
                  'Catalan_convergence': ratios,
                  'tree_classes': 14, 'one_labelling_per_class': True,
                  'scope': 'Normalization check; not a re-certification of implicit-function entries.'}, indent=2))