"""Independent numerical and stored-interval checks for numberdb.org/T611.

    sage -python checks.py
    sage -python checks.py --stored

Requires numberdb and mpmath. --stored requires a key while T611 is private.
The generator uses Arb Li_2; this file uses mpmath quadrature, without Li_2.
Antiprisms: integrate the defining log-sine integral independently.
Lobell: use the geometrically derived arccosh integral in Mednykh--Vesnin,
Lobell manifolds revised (2007), Theorem 2, divided by eight for a polyhedron.
This tests the volume formula, not only its special-function evaluation.

A sign discrepancy is recorded deliberately: equation (2) in that paper has
+Lambda(2*theta-pi/2), whereas Inoue Theorem 4.1 has minus. The former sign
contradicts the same paper's independently derived Theorem 2 integral and
known dodecahedral volume. We use Inoue's formula and verify the discrepancy.

Numerical quadrature is a cross-check, not the proof of 100 certified digits.
Controls run first: integral 1, log-sine at pi/4 vs Catalan, and independent
arccosh integral at alpha=pi/4 vs Catalan/2. All comparisons require finite
answers and stated tolerances. The same checks apply to every database row.
"""
import json
import os
import sys
from decimal import Decimal
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
import mpmath as mp
from generate import RightAngledVolumes, volume


def log_sine(x):
    # Subtract the endpoint logarithm. The ratio has the continuous value 1.
    def f(u):
        if not u:
            return mp.mpf(0)
        return mp.log(mp.sin(x*u)/(x*u))
    return x*(1-mp.log(2*x)) - x*mp.quad(f,[0,1])


def schlafli_integral(alpha):
    # Split at the logarithmic singularity t=pi/4.
    # Its explicit primitive prevents lost digits near the singular endpoint.
    c=mp.pi/4
    b=mp.pi/3
    endpoint_tolerance=100*mp.eps
    def h(t):
        u=t-c
        ratio=1 if not u else mp.sin(2*u)/(2*u)
        radicand=mp.cos(t)**2-mp.cos(2*t)**2
        # The integrand is real on [alpha,b]. mpmath can have an endpoint
        # imaginary roundoff of order sqrt(eps); take the real part there.
        assert radicand >= -endpoint_tolerance
        return mp.log(mp.cos(t)+mp.sqrt(max(mp.mpf(0),radicand)))-mp.log(ratio)
    def primitive(u):
        return mp.mpf(0) if not u else u*(mp.log(2*abs(u))-1)
    return mp.quad(h,[alpha,c,b]) - primitive(b-c)+primitive(alpha-c)


def independent(family,n):
    if family=='antiprism':
        return 2*n*(log_sine(mp.pi/4+mp.pi/(2*n))
                    +log_sine(mp.pi/4-mp.pi/(2*n)))
    return n*schlafli_integral(mp.pi/n)


def decimal_interval(text):
    value=Decimal(text)
    ulp=Decimal(1).scaleb(value.as_tuple().exponent)
    return mp.mpf(str(value)),mp.mpf(str(ulp))


def rational_decimal(text):
    value=Decimal(text)
    sign,digits,exponent=value.as_tuple()
    integer=int(''.join(map(str,digits)))*(-1 if sign else 1)
    return QQ(integer)*QQ(10)**exponent


def check_decimal(text,ball,independent_value):
    centre,ulp=decimal_interval(text)
    assert mp.isfinite(independent_value)
    error=abs(centre-independent_value)
    assert error < ulp, (text,mp.nstr(error/ulp,10))
    # Check exact decimal endpoints contain the certified computation.
    exact_centre=rational_decimal(text)
    exact_ulp=QQ(10)**Decimal(text).as_tuple().exponent
    real=RealBallField(600)
    assert real(exact_centre-exact_ulp) < ball < real(exact_centre+exact_ulp)
    return error/ulp


def stored_rows():
    document=numberdb.table('T611')
    block=document['Numbers']
    if isinstance(block,list):
        return [(str(r['params']['family']),int(r['params']['n']),r['number']) for r in block]
    rows=[]
    for family,entries in block.items():
        entries=entries.get('numbers',entries)
        for n,entry in entries.items():
            if n=='param-latex':continue
            rows.append((family,int(n),entry['number'] if isinstance(entry,dict) else entry))
    return rows


def run(use_stored=False):
    mp.mp.dps=135
    tol=mp.mpf('1e-120')
    controls={
        'quadrature_constant':abs(mp.quad(lambda x:mp.mpf(1),[0,1])-1),
        'log_sine_Catalan':abs(log_sine(mp.pi/4)-mp.catalan/2),
        'Schlafli_ideal_control':abs(schlafli_integral(mp.pi/4)-mp.catalan/2)}
    assert all(mp.isfinite(x) and x < tol for x in controls.values()),controls
    generator=RightAngledVolumes()
    parameters=[(p['family'],p['n']) for p in generator.enumerate()]
    assert parameters==([('antiprism',n) for n in range(3,31)]
                         +[('lobell',n) for n in range(5,31)])
    balls={(family,n):volume(family,n) for family,n in parameters}
    texts={(family,n):to_text(b, digits=100) for (family,n),b in balls.items()}
    if use_stored:
        rows=stored_rows()
        assert [(f,n) for f,n,t in rows]==parameters
        texts={(f,n):t for f,n,t in rows}
    output=[]
    max_ulp=mp.mpf(0)
    min_digits=1000.0
    for family,n in parameters:
        ball=balls[family,n]
        expected=independent(family,n)
        text=texts[family,n]
        ulp_error=check_decimal(text,ball,expected)
        max_ulp=max(max_ulp,ulp_error)
        relative_radius=float(ball.rad())/abs(float(ball.mid()))
        import math
        min_digits=min(min_digits,-math.log10(relative_radius))
        # Independent quadrature at a second precision detects instability.
        with mp.workdps(115):
            second=independent(family,n)
        assert abs(expected-second) < mp.mpf('1e-105')*abs(expected)
        output.append({'family':family,'n':n,'number':text,
                       'quadrature':mp.nstr(expected,120),
                       'error_in_last_place':mp.nstr(ulp_error,12)})
    assert abs(independent('antiprism',3)-4*mp.catalan) < tol
    assert abs(independent('antiprism',4)-independent('lobell',6)) < tol
    # Independently published 100 decimal digits of L_5 (OEIS A321209).
    known=mp.mpf('4.306207600730808652919837159842209449241185235412934114224949117819904343308238335156704866571730974')
    assert abs(independent('lobell',5)-known) < mp.mpf('1e-92')
    # Demonstrate the printed sign error in Mednykh--Vesnin equation (2).
    n=5
    theta=mp.asin(1/(2*mp.cos(mp.pi/n)))
    incorrect=mp.mpf(n)/2*(2*log_sine(theta)+log_sine(theta+mp.pi/n)
                           +log_sine(theta-mp.pi/n)-log_sine(mp.pi/2-2*theta))
    correct=independent('lobell',5)
    assert abs(correct-incorrect)>1
    return {'table':'T611','stored':use_stored,'entries':len(output),
            'controls':{k:mp.nstr(v,8) for k,v in controls.items()},
            'min_relative_radius_decimal_digits':min_digits,
            'max_error_in_last_place':mp.nstr(max_ulp,12),
            'quadrature_precisions':[135,115],
            'second_precision_relative_tolerance':'1e-105',
            'source_sign_check':{'Inoue_and_integral':mp.nstr(correct,30),
                                 'Mednykh_Vesnin_equation_2_printed_sign':mp.nstr(incorrect,30)},
            'values':output,
            'skill_manifest':'agents/runs/20261029T183732Z-build-skills/manifest.json',
            'skill_sha256_local_and_live':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
            'mpmath_version':mp.__version__}


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