"""Independent checks for T530; run with sage -python checks.py [--stored].

The resolvent method is derived in derivation.md. It shares neither the
K,Q,P,L recurrences nor the parity simplification with generate.py. Exact
Airy-series division provides a third coefficient computation. The repository
helpers in check.py supply exactness, agreement, ring and size checks.
--stored additionally fetches the draft and checks actual written intervals;
set NUMBERDB_API_KEY or NUMBERDB_KEY_FROM_STDIN=1 for a private draft.
"""
import os
import sys
import json
import hashlib
from pathlib import Path
from math import factorial
from fractions import Fraction
from decimal import Decimal
import numberdb.sage as numberdb
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
import mpmath as mp
from check import exactness, agrees_with, names_its_rings, measure
from one_value import moment
from generate import PositiveBrownianMoments, coefficients, rational_factor, MAX_ORDER


def derivative(poly):
    return [i*poly[i] for i in range(1, len(poly))] or [QQ(0)]


def add(a, b):
    return [(a[i] if i < len(a) else QQ(0))
            + (b[i] if i < len(b) else QQ(0))
            for i in range(max(len(a),len(b)))]


def scale(poly, c):
    return [c*x for x in poly]


def resolvent(last):
    """V_n(x)=P_n(x)+exp(-x)Q_n(x), for x>=0 and Exp(1/2) time."""
    polynomial = [QQ(1)]
    exponential = [QQ(0)]
    constants = [QQ(1)]
    for n in range(1,last+1):
        rhs = [QQ(0)] + scale(polynomial, QQ(2*n))
        solution = [QQ(0)]
        term = rhs
        while any(term):
            solution = add(solution, term)
            term = derivative(derivative(term))
        assert add(solution, scale(derivative(derivative(solution)),QQ(-1))) == rhs
        rhs_exp = [QQ(0)] + scale(exponential, QQ(2*n))
        derivative_solution = [QQ(0)]
        term = scale(rhs_exp, QQ(1)/2)
        while any(term):
            derivative_solution = add(derivative_solution, term)
            term = scale(derivative(term), QQ(1)/2)
        # Continuity of V_n and V_n' at zero matches C_n exp(x) at x<0.
        constant = (derivative(solution)[0] + derivative_solution[0] - solution[0])/QQ(2)
        exp_solution = [constant] + [x/QQ(i+1) for i,x in enumerate(derivative_solution)]
        check_exp = add(scale(derivative(exp_solution),QQ(2)),
                        scale(derivative(derivative(exp_solution)),QQ(-1)))
        # Trailing zero coefficients do not change a polynomial.
        def normalized(a):
            a = a[:]
            while len(a)>1 and a[-1]==0: a.pop()
            return a
        assert normalized(check_exp) == normalized(rhs_exp)
        assert derivative(solution)[0]+derivative(exp_solution)[0]-exp_solution[0] == solution[0]+exp_solution[0]
        polynomial, exponential = solution, exp_solution
        constants.append(polynomial[0]+exponential[0])
        if n == 1:
            assert normalized(polynomial) == [QQ(0),QQ(2)]
            assert normalized(exponential) == [QQ(1)]
    assert constants[:3] == [QQ(1),QQ(1),QQ(17)/2]
    return constants


def airy_quotient(last):
    alpha = [QQ(factorial(6*n))/(QQ(576)**n*factorial(3*n)*factorial(2*n))
             for n in range(last+1)]
    beta = [QQ(1)]
    for n in range(1,last+1):
        beta.append(alpha[n]+QQ(3*(2*n-1))/4*beta[n-1])
    denominator = [QQ(2)*alpha[n]/QQ(1-6*n) for n in range(last+1)]
    result = []
    for n in range(last+1):
        result.append((alpha[n]+beta[n]-sum(denominator[j]*result[n-j]
                       for j in range(1,n+1)))/denominator[0])
    assert result[:3] == [QQ(1),QQ(1)/2,QQ(17)/16]
    return result


def exact_decimal(text):
    f = Fraction(text)
    return QQ(f.numerator)/QQ(f.denominator)


def written_bounds(text):
    centre = exact_decimal(text)
    if '.' not in text and 'e' not in text.lower():
        return centre, centre
    unit = QQ(10)**int(Decimal(text).as_tuple().exponent)
    return centre-unit, centre+unit


def check_written(text, expected, exact=None):
    lo,hi = written_bounds(text)
    if exact is not None:
        assert lo==hi==exact
    else:
        assert lo<hi and expected.is_finite()
        assert expected.accuracy() > 700
        assert expected > lo and expected < hi


def gaussian_checks():
    mp.mp.dps = 85
    control = mp.quad(lambda t: mp.sin(t)**2*mp.cos(t)**2, [0,mp.pi/2])
    assert abs(control-mp.pi/16) < mp.mpf('1e-80')
    mean = mp.quad(lambda t: mp.sqrt(t/(2*mp.pi)),[0,1])
    second = 2/(3*mp.pi)*mp.quad(
        lambda t: mp.sin(t)**2*(mp.cos(t)+mp.sin(t)*(mp.pi/2+t))*mp.cos(t),
        [0,mp.pi/2])
    assert abs(mean-mp.sqrt(2/mp.pi)/3)<mp.mpf('1e-80')
    assert abs(second-mp.mpf(17)/96)<mp.mpf('1e-80')
    return {'working_decimal_digits':85,'nonzero_control':'integral sin^2(theta) cos^2(theta) = pi/16',
            'control_error':mp.nstr(abs(control-mp.pi/16),6),
            'mean':mp.nstr(mean,80),'second_moment':mp.nstr(second,80),
            'method':'Gaussian marginal mean and joint positive-part product integral; agreement check, not certification'}


def run(stored=False):
    base = Path(__file__).resolve().parent
    gen = PositiveBrownianMoments()
    params = list(gen.enumerate())
    assert params == [{'n':n} for n in range(1,13)]
    assert not names_its_rings(str(base/'generate.py'))
    values = {p['n']:gen.value(p,100) for p in params}
    assert not exactness(values)
    gaussian = gaussian_checks()  # nonzero control precedes this measurement
    fk = resolvent(MAX_ORDER)
    ell = coefficients(MAX_ORDER)[3]
    assert not agrees_with(dict(enumerate(ell)),lambda n:airy_quotient(MAX_ORDER)[n])
    assert not agrees_with(dict(enumerate(ell)),lambda n:fk[n]/(QQ(2)**n*factorial(n)))
    assert rational_factor(0)==1
    for n,c in {1:QQ(1)/3,2:QQ(17)/96,3:QQ(251)/1260,4:QQ(6989)/40960}.items():
        assert rational_factor(n)==c
    high = RealBallField(768)
    independent = {n:high(fk[n])/(high(2)**(QQ(3*n)/2)*(high(QQ(3*n)/2+1)).gamma())
                   for n in range(MAX_ORDER+1)}
    # Perman and Wellner, printed Table 1, positive-motion column. Odd
    # coefficients multiply 1/sqrt(2*pi), not sqrt(2/pi).
    printed = {1:QQ(2)/3, 2:QQ(17)/96, 3:QQ(251)/630,
               4:QQ(6989)/40960, 5:QQ(188267)/360360,
               6:QQ(37235311)/132120576, 7:QQ(1451280043)/1396755360,
               8:QQ(522258818027)/797253304320, 9:QQ(8096107769)/2930256000,
               10:QQ(1011667945773427)/515911471595520}
    for n,c in printed.items():
        assert rational_factor(n)==(c/2 if n%2 else c)
    rows = []
    accuracy = []
    for n,v in values.items():
        exact = fk[n]/(QQ(2)**(3*n//2)*factorial(3*n//2)) if n%2==0 else None
        if n%2:
            assert v.is_finite() and v.accuracy()>=500
            accuracy.append(int(v.accuracy()))
            assert v.overlaps(independent[n])
        else:
            assert v==exact and (v in independent[n])
        program_value = moment(n)
        if n%2:
            assert program_value.is_finite() and program_value.accuracy()>490
            assert program_value.overlaps(independent[n])
        else:
            assert program_value==exact
        text = to_text(v,digits=100)
        check_written(text,independent[n],exact)
        rows.append({'params':{'n':str(n)},'number':text})
    for n in range(1,MAX_ORDER):
        assert independent[n-1]*independent[n+1]-independent[n]**2 > 0
    all_recurrences = True
    k,q,p,l = coefficients(MAX_ORDER)
    for n in range(1,MAX_ORDER+1):
        assert k[n]==QQ(3*n-4)/4*k[n-1]+sum(k[j]*k[n-j] for j in range(1,n))
        assert q[n]==QQ(3*n-2)/2*q[n-1]+2*sum(k[j]*q[n-j] for j in range(1,n+1))
        assert p[n]==sum(k[j]*p[n-j] for j in range(1,n+1))
        assert 2*l[n]==p[n]+sum(p[j]*q[n-j] for j in range(n+1))
    stored_report = None
    if stored:
        doc = numberdb.table('T530')
        records = _stored_entries(doc)
        assert len(records)==12
        assert list(records)==[str(n) for n in range(1,13)]
        balls = {0: high(1)}
        for identity,text in records.items():
            n=int(identity)
            exact = fk[n]/(QQ(2)**(3*n//2)*factorial(3*n//2)) if n%2==0 else None
            check_written(text,independent[n],exact)
            lo,hi=written_bounds(text)
            balls[n]=high((lo+hi)/2).add_error(high((hi-lo)/2))
            # Invert the moment formula on the stored interval, then test
            # the independently computed coefficient and convolution.
            recovered = balls[n]*(high(QQ(3*n)/2+1)).gamma()*high(2)**(QQ(n)/2)/factorial(n)
            assert (ell[n] in recovered)
            assert (airy_quotient(MAX_ORDER)[n] in recovered)
            assert recovered.accuracy()>320
        for n in range(1,MAX_ORDER):
            assert balls[n-1]*balls[n+1]-balls[n]**2 > 0
        stored_report={'entries':12,'all_independent_enclosures_contained':True,
                       'inverted_moment_formula_and_convolution':True,'strict_log_convexity':True}
    measured=measure({str(n):to_text(v,100) for n,v in values.items()})
    evidence={'table':'T530','orders':[1,12],'entries':12,'exact_even_moments':6,
              'certified_odd_moments':6,'stored_digits':100,'working_bits':numberdb.bits(100,losing=180),
              'minimum_relative_accuracy_bits':min(accuracy),'per_odd_order_accuracy_bits':accuracy,
              'independent_working_bits':768,'resolvent_coefficients':[str(x) for x in fk],
              'rational_factors':[str(rational_factor(n)) for n in range(1,13)],
              'airy_quotient_all_orders':True,'feynman_kac_all_orders':True,
              'all_stated_recurrences':all_recurrences,'janson_table8_orders_1_to_4':True,
              'perman_wellner_table1_orders_1_to_10':True,'program_checked_all_orders':True,
              'written_interval_containment':True,'strict_log_convexity_orders_1_to_11':True,
              'gaussian_checks':gaussian,'measurement':measured,'stored_readback':stored_report,
              'instruction_manifest':'agents/runs/20261026T220819Z-build-skills/manifest.json',
              'instruction_sha256':'b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18',
              'generator_sha256':hashlib.sha256((base/'generate.py').read_bytes()).hexdigest()}
    from sage.version import version as sage_version
    import importlib.metadata
    evidence['versions']={'sage':sage_version,'mpmath':mp.__version__,
                          'numberdb':importlib.metadata.version('numberdb')}
    print('EVIDENCE_JSON='+json.dumps(evidence,sort_keys=True))
    print('ROWS_JSON='+json.dumps(rows))
    return evidence


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    run(stored='--stored' in sys.argv)