"""Independent checks for Hessian Mahler measures -- numberdb.org/T649.

Run: sage -python checks.py (with generate.py and values.json alongside).
Numerical checks use the original cubic, not the quadratic integrand used
by generate.py. At each x=e^(i theta), solve the cubic in y by Cardano's
formula and integrate sum(log(max(1,abs(y)))) over 0<=theta<=pi/3.
Symmetry gives the factor 3/pi. Split at all projected torus zeros and at
Cardano discriminant zeros; midpoint splits aid convergence. This check
is also floating quadrature, not certification.
"""
import json
import time
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
import mpmath as mp
from generate import ITEMS,mp_argument,exact_argument


def cubic_jensen(item,working=125):
    mp.mp._tanh_sinh.clear()
    with mp.workdps(working):
        t=mp_argument(item)
        k=mp.sign(t)*mp.root(abs(t),3)
        omega=mp.exp(2*mp.pi*mp.j/3)
        def integrand(theta):
            x=mp.exp(mp.j*theta)
            p=-k*x
            q=x**3+1
            if t==0:
                return max(mp.mpf(0),mp.log(abs(q)))
            if t==27:
                ys=[-x-1,-omega*x-omega**2,-omega**2*x-omega]
                return sum(max(mp.mpf(0),mp.log(abs(y))) for y in ys)
            root=mp.sqrt(q*q/4+p**3/27)
            choices=(-q/2+root,-q/2-root)
            u=max(choices,key=abs)**(mp.mpf(1)/3)
            v=-p/(3*u)
            ys=[u+v,omega*u+omega**2*v,omega**2*u+omega*v]
            return sum(max(mp.mpf(0),mp.log(abs(y))) for y in ys)
        end=mp.pi/3
        points=[mp.mpf(0),end]
        if -1<t<27:
            c=mp.acos((k-1)/2)
            points.extend([c/3,min(2*c/3,(2*mp.pi-2*c)/3)])
        if 0<t<27:
            points.append(2*mp.acos(mp.sqrt(t/27))/3)
        points=sorted(set(points))
        # Remove coincident points within working-precision rounding.
        unique=[points[0]]
        for p in points[1:]:
            if p-unique[-1]>mp.eps*100:
                unique.append(p)
        split=[unique[0]]
        for a,b in zip(unique,unique[1:]):split.extend([(a+b)/2,b])
        tolerance=mp.power(10,-working+15)
        def adaptive(a,b,depth=0):
            value,error=mp.quad(integrand,[a,b],error=True,maxdegree=8)
            if error<tolerance*(b-a)/end:
                return value,error
            if depth>=24:raise ArithmeticError('Cubic check convergence: '+item['key'])
            m=(a+b)/2
            left,el=adaptive(a,m,depth+1);right,er=adaptive(m,b,depth+1)
            return left+right,el+er
        pieces=[adaptive(a,b) for a,b in zip(split,split[1:])]
        result=sum(v for v,e in pieces);error=sum(e for v,e in pieces)
        assert mp.isfinite(result)
        return mp.nstr(result*3/mp.pi,working),mp.nstr(error*3/mp.pi,8)


def hypergeometric(item,working=125):
    with mp.workdps(working):
        t=mp_argument(item)
        assert t<=-1 or t>=27
        if abs(t)>27:
            # Constant-term hypergeometric series. Its nth term is
            # (3n)! / (3n*(n!)^3*t^n). The absolute ratio is <27/abs(t),
            # so the remaining absolute tail is bounded geometrically.
            # Explicit iteration avoids mpmath's expensive acceleration
            # near the unit argument. Rounding remains non-certified.
            term=2/t
            total=term
            n=1
            ratio_bound=27/abs(t)
            target=mp.power(10,-working+5)
            while True:
                term*=mp.mpf(3*n*(3*n+1)*(3*n+2))/((n+1)**3*t)
                n+=1
                total+=term
                if abs(term)*ratio_bound/(1-ratio_bound)<target:break
                if n>1000000:raise ArithmeticError('Series did not converge')
            val=mp.log(abs(t))/3-total
        else:
            val=mp.log(abs(t))/3-mp.re(2/t*mp.hyper([mp.mpf(4)/3,mp.mpf(5)/3,1,1],[2,2,2],27/t))
        return mp.nstr(val,working)


def agreement(a,b):
    with mp.workdps(150):
        delta=abs(mp.mpf(a)-mp.mpf(b))/max(1,abs(mp.mpf(b)))
        return float(-mp.log10(delta)) if delta else 999.0


def symbolic_checks():
    # c4^3/discriminant from the shared Weierstrass model.
    R=PolynomialRing(QQ,'t');t=R.gen()
    a2=-27*t**2;a4=216*t**3*(t-27);a6=-432*t**4*(t-27)**2
    b2=4*a2;b4=2*a4;b6=4*a6;b8=4*a2*a6-a4*a4
    c4=b2*b2-24*b4
    disc=-b2*b2*b8-8*b4**3-27*b6**2+9*b2*b4*b6
    assert c4**3*(t-27)**3==disc*t*(t+216)**3
    # Torus monomial substitution, after clearing negative powers:
    # (x^2*y)^3 P(y/x^2,1/(x*y)) = Q(x^3,y^3).
    R3=PolynomialRing(QQ,('x','y','k'));x,y,k=R3.gens()
    assert y**6+x**3+x**6*y**3-k*x**3*y**3 == y**6+(x**6-k*x**3)*y**3+x**3


def run(values):
    symbolic_checks()
    # A known, nonzero control for the independent integrator runs first.
    with mp.workdps(130):
        L=(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
        control=mp.nstr(3*mp.sqrt(3)*L/(4*mp.pi),125)
    zero=next(i for i in ITEMS if i['key']=='0')
    c,err=cubic_jensen(zero)
    assert agreement(c,control)>115
    print(json.dumps({'control':'N(0) from Hurwitz zeta Dirichlet L-value','digits':agreement(c,control),'value':c}),flush=True)
    for item in ITEMS:
        key=item['key'];start=time.monotonic()
        v,err=cubic_jensen(item)
        digs=agreement(v,values[key]['high'])
        assert digs>110,(key,digs,v,values[key]['high'])
        row={'t':key,'cubic_jensen':v,'quadrature_estimate':err,'agreement_digits':digs}
        if exact_argument(item)<=-1 or exact_argument(item)>27:
            h=hypergeometric(item)
            hd=agreement(h,values[key]['high'])
            assert hd>110,(key,hd)
            row.update(hypergeometric=h,hypergeometric_agreement_digits=hd)
        if key=='27':
            with mp.workdps(130):target=mp.nstr(3*mp.mpf(control),125)
            assert agreement(v,target)>115
        row['seconds']=round(time.monotonic()-start,2)
        print(json.dumps(row),flush=True)


if __name__=='__main__':
    values=json.loads(Path(__file__).with_name('values.json').read_text())
    run(values)