#!/usr/bin/env python3
# Author: David Niedbala Giraudin (ORCID 0009-0009-1526-1178), independent researcher, France.
# Ancillary file of the note 'The Schur-Siegel-Smyth trace problem: the lower bound 1.8022246'.
# Written with the assistance of an AI system (Claude, Anthropic).
"""Independent rigorous checker for lower-bound certificates of the Schur-Siegel-Smyth trace constant.

Lemma (Orloski-Talebizadeh Sardari-Smith; Prop. 2 of arXiv:2609.29298; Lemma 2.1 of Tomar):
C_86 >= lam if  F(x) = x - lam - sum_Q lam_Q log|Q(x)| - 2 U_g(x) + I(g)/s >= 0  for all x>0, x not a root,
where g >= 0 is a bounded density with compact support and mass s>0, U_g(x)=int log|x-y| g(y)dy, I(g)=int U_g g.

Certificate format (Tomar, CERT_FORMAT.md): A (integer polys, highest degree first), lamQ (rationals),
nodes e_0<...<e_N, rho (N-1 rationals; rho[k-1] = weight of the unit-mass hat on [e_{k-1},e_{k+1}]),
lambda, x_lo, x_hi.

Method (written from the lemma only):
  g is continuous piecewise linear, g'' = sum_n C_n delta_{e_n} with exact rational C_n, hence
  U_g(x) = sum_n C_n G(x-e_n),  G(t) = t^2 log|t|/2 - 3t^2/4   (G''=log|t|),
  I(g)   = sum_{n,m} C_n C_m Psi(e_n-e_m), Psi(t) = t^4 log|t|/24 - 25 t^4/288  (Psi''''=log|t|).
  All arithmetic in arb balls (python-flint). On [x_lo, x_hi] the line is cut at all nodes, at isolating
  intervals of the real roots, and into blocks; on each block the contribution of distant nodes is a Taylor
  polynomial with an explicit remainder bound; on each cell J (no node inside) F is bounded below by
  F(m)+F'(m)(x-m) - max(0,-inf_J F''_rest)(w/2)^2/2 plus chord/tangent bounds for the singular endpoint terms.
  Tails (0,x_lo] and [x_hi,inf) by monotonicity bounds.
"""
import sys, json, time
from fractions import Fraction
import flint
from flint import arb, fmpq, fmpz_poly, fmpq_poly, ctx

PREC = 192
ctx.prec = PREC

def Fr(s): return Fraction(s)
def toarb(q): return arb(fmpq(q.numerator, q.denominator))

# ---------- exact data ----------
def load(path):
    c = json.load(open(path))
    A = [[int(a) for a in Q] for Q in c['A']]
    lamQ = [Fr(x) for x in c['lamQ']]
    e = [Fr(x) for x in c['nodes']]
    rho = [Fr(x) for x in c['rho']]
    lam = Fr(c['lambda']); xlo = Fr(c['x_lo']); xhi = Fr(c['x_hi'])
    assert len(rho) == len(e) - 2
    assert all(e[i] < e[i+1] for i in range(len(e)-1)), "nodes not increasing"
    assert all(r >= 0 for r in rho) and all(l >= 0 for l in lamQ), "negative weight"
    return A, lamQ, e, rho, lam, xlo, xhi

def jumps(e, rho):
    """C_n = jump of g' at e_n (exact rationals), n=0..N; g(e_k)=rho_k*2/(e_{k+1}-e_{k-1})."""
    N = len(e) - 1
    gv = [Fraction(0)] * (N + 1)
    for k in range(1, N):
        gv[k] = rho[k-1] * 2 / (e[k+1] - e[k-1])
    slope = [(gv[k+1] - gv[k]) / (e[k+1] - e[k]) for k in range(N)]
    C = [Fraction(0)] * (N + 1)
    for n in range(N + 1):
        left = slope[n-1] if n >= 1 else Fraction(0)
        right = slope[n] if n < N else Fraction(0)
        C[n] = right - left
    return C, gv

# ---------- special functions in arb ----------
def _small(t):
    # ball t containing 0 (e.g. difference of two enclosures of the same rational): return eps=sup|t| (< 1)
    eps = abs(t).upper()
    assert eps < 1e-30
    return eps
def G(t):        # t^2 log|t|/2 - 3t^2/4 ;  |G(t)| <= eps^2(|log eps|/2+3/4) for |t|<=eps<1
    if t.contains(0):
        eps = _small(t)
        if eps == 0: return arb(0)
        return arb(0).union(eps*eps*(-eps.log()/2 + arb(3)/4)).union(-eps*eps*(-eps.log()/2 + arb(3)/4))
    return t*t*(abs(t).log()/2 - arb(3)/4)
def G1(t):       # t log|t| - t ; |G1| <= eps(|log eps|+1)
    if t.contains(0):
        eps = _small(t)
        if eps == 0: return arb(0)
        b = eps*(-eps.log() + 1)
        return b.union(-b)
    return t*(abs(t).log() - 1)
def Psi(t):
    if t.contains(0): raise ValueError("Psi at 0 should not occur")
    t4 = t**4
    return t4*abs(t).log()/24 - arb(25)/288*t4

# ---------- energy ----------
def energy(e, C):
    """I(g) = sum_{n,m} C_n C_m Psi(e_n - e_m). Polynomial part exact: -25/288 * 6 M_2^2 (since M_0=M_1=0)."""
    M = [sum(Cn * en**j for Cn, en in zip(C, e)) for j in range(5)]
    assert M[0] == 0 and M[1] == 0, "g'' moments M0,M1 must vanish"
    poly_part = Fraction(-25, 288) * (2*M[4]*M[0] - 8*M[3]*M[1] + 6*M[2]**2)
    Ca = [toarb(c) for c in C]; ea = [toarb(x) for x in e]
    nz = [i for i, c in enumerate(C) if c != 0]
    S = arb(0)
    for a in range(len(nz)):
        i = nz[a]; ci = Ca[i]; ei = ea[i]; acc = arb(0)
        for b in range(a+1, len(nz)):
            j = nz[b]; t = ei - ea[j]; t2 = t*t
            acc += Ca[j]*t2*t2*abs(t).log()
        S += ci*acc
    log_part = 2*S/24
    return toarb(poly_part) + log_part

# ---------- roots ----------
def real_roots_isolated(Q, width=Fraction(1, 2**90)):
    """All roots real (checked by Sturm count) and isolated in rational intervals of tiny width."""
    p = fmpq_poly(list(reversed(Q)))
    d = p.degree()
    # Sturm sequence
    seq = [p, p.derivative()]
    while seq[-1].degree() > 0:
        r = -(seq[-2] % seq[-1])
        if r == 0: raise ValueError("repeated roots")
        seq.append(r)
    def sign_changes(x):
        vals = [q(x) for q in seq]
        vals = [v for v in vals if v != 0]
        return sum(1 for a, b in zip(vals, vals[1:]) if (a > 0) != (b > 0))
    def sc_inf(sgn):
        vals = []
        for q in seq:
            lc = q.coeffs()[-1] if q.degree() >= 0 else 0
            vals.append(lc * (sgn ** q.degree()))
        vals = [v for v in vals if v != 0]
        return sum(1 for a, b in zip(vals, vals[1:]) if (a > 0) != (b > 0))
    nreal = sc_inf(-1) - sc_inf(1)
    assert nreal == d, f"polynomial {Q} is not real-rooted"
    # bound roots
    co = [abs(Fraction(int(c))) for c in Q]
    B = 1 + max(co[1:]) / co[0] if len(co) > 1 else Fraction(1)
    out = []
    def rec(a, b, cnt):
        if cnt == 0: return
        if cnt == 1 and b - a <= width:
            out.append((a, b)); return
        m = (a + b) / 2
        fm = fmpq(m.numerator, m.denominator)
        if p(fm) == 0:
            out.append((m, m))
            eps = width / 8
            l = m - eps; r = m + eps
            rec(a, l, sign_changes(fmpq(a.numerator, a.denominator)) - sign_changes(fmpq(l.numerator, l.denominator)))
            rec(r, b, sign_changes(fmpq(r.numerator, r.denominator)) - sign_changes(fmpq(b.numerator, b.denominator)))
            return
        ca = sign_changes(fmpq(a.numerator, a.denominator)); cmid = sign_changes(fm); cb = sign_changes(fmpq(b.numerator, b.denominator))
        rec(a, m, ca - cmid); rec(m, b, cmid - cb)
    a0 = -B - 1; b0 = B + 1
    rec(a0, b0, sign_changes(fmpq(a0.numerator, a0.denominator)) - sign_changes(fmpq(b0.numerator, b0.denominator)))
    out.sort()
    assert len(out) == d
    return out

# ---------- main pointwise check ----------
def fl(q): return float(q)

class Checker:
    def __init__(self, path, P=24, L=Fraction(1, 50)):
        t0 = time.time()
        self.A, self.lamQ, self.e, self.rho, self.lam, self.xlo, self.xhi = load(path)
        self.C, self.gv = jumps(self.e, self.rho)
        self.s = sum(self.rho)
        self.ea = [toarb(x) for x in self.e]; self.Ca = [toarb(c) for c in self.C]
        self.ef = [fl(x) for x in self.e]
        self.P = P; self.L = L
        # roots of polynomials with positive multiplier
        self.R = []   # (lamQ arb, lo Fraction, hi Fraction, arb ball)
        self.leadterm = arb(0)
        for Q, l in zip(self.A, self.lamQ):
            if l == 0: continue
            la = toarb(l)
            self.leadterm += la * arb(abs(Q[0])).log()
            for lo, hi in real_roots_isolated(Q):
                ball = toarb(lo).union(toarb(hi))
                assert ball.contains(toarb(lo)) and ball.contains(toarb(hi))
                self.R.append((la, lo, hi, ball))
        self.sumlamdeg = sum(l * (len(Q) - 1) for Q, l in zip(self.A, self.lamQ))
        print(f"loaded: {len(self.A)} polys, N={len(self.e)-1} cells, s={fl(self.s):.12f}, lambda={self.lam} = {fl(self.lam):.12f}  ({time.time()-t0:.1f}s)", flush=True)

    def setup_energy(self):
        t0 = time.time()
        self.I = energy(self.e, self.C)
        self.Is = self.I / toarb(self.s)
        print(f"I(g) = {self.I}   I(g)/s = {self.Is}  ({time.time()-t0:.0f}s)", flush=True)
        self.const = -toarb(self.lam) + self.Is - self.leadterm   # Lin(x) = x + const

    def Ug_point(self, x):  # rigorous U_g at a point (all nodes), for tails/controls
        S = arb(0)
        for c, en in zip(self.Ca, self.ea):
            if not c.is_zero(): S += c * G(x - en)
        return S

    def tails(self):
        ok = True
        # left: (0, x_lo]; need every root with positive multiplier to be either 0 or > x_lo
        xl = toarb(self.xlo)
        bound = self.const  # x >= 0 dropped (x - ...): x>0 adds positive
        for la, lo, hi, ball in self.R:
            if hi == 0 and lo == 0:
                bound += -la * xl.log()          # -lam log x >= -lam log x_lo on (0,x_lo]
            else:
                assert lo > self.xlo, "root in (0,x_lo]"
                bound += -la * ball.log()        # -lam log(r-x) >= -lam log r
        U0 = arb(0)
        nzk = [k for k in range(len(self.gv)) if self.gv[k] != 0]
        supp_lo = self.e[nzk[0] - 1]; supp_hi = self.e[nzk[-1] + 1]
        self.supp = (supp_lo, supp_hi)
        assert supp_lo > 0 and self.xlo <= supp_lo, "x_lo must not exceed the start of supp g"
        U0 = self.Ug_point(arb(0))              # U_g(x) <= U_g(0) for x in (0, supp_lo]
        bound += -2 * U0
        okL = bound > 0
        print(f"left tail (0,{self.xlo}]: {okL}  bound {bound}", flush=True)
        # right: [x_hi, inf): F >= x - beta log x + const, beta = sum lam_Q deg Q + 2s, increasing for x>beta
        beta = self.sumlamdeg + 2 * self.s
        assert self.xhi >= beta and self.xhi > self.supp[1]
        for la, lo, hi, ball in self.R: assert hi < self.xhi and lo >= 0
        xh = toarb(self.xhi)
        val = xh - toarb(beta) * xh.log() + self.const
        okR = val > 0
        print(f"right tail [{self.xhi},inf): {okR}  bound {val}  beta={fl(beta):.6f}", flush=True)
        return okL and okR

    def block_taylor(self, c, r, far):
        """coefficients a_k of sum_far C_n G(x-e_n) about c, and remainder bound for |x-c|<=r."""
        P = self.P; ca = toarb(c)
        a = [arb(0) for _ in range(P + 1)]
        rem = arb(0)
        ra = toarb(r)
        fact = [1]
        for k in range(1, P + 2): fact.append(fact[-1] * k)
        for n in far:
            Cn = self.Ca[n]
            if Cn.is_zero(): continue
            t = ca - self.ea[n]
            a[0] += Cn * G(t); a[1] += Cn * G1(t); a[2] += Cn * abs(t).log() / 2
            inv = 1 / t; pw = inv  # 1/t^{k-2}
            for k in range(3, P + 1):
                a[k] += Cn * ((-1) ** (k - 3)) * arb(fact[k - 3]) / arb(fact[k]) * pw
                pw *= inv
            d = abs(t); q = ra / d
            rem += abs(Cn) * d * d * q ** (P + 1) / arb((P + 1) * P * (P - 1)) / (1 - q)
        return a, rem

    def run(self, maxdepth=40, verbose_fail=20):
        t0 = time.time()
        okT = self.tails()
        # breakpoints
        bps = set(x for x in self.e if self.xlo < x < self.xhi)
        bps.add(self.xlo); bps.add(self.xhi)
        for la, lo, hi, ball in self.R:
            for z in (lo, hi):
                if self.xlo < z < self.xhi: bps.add(z)
        # block grid: fine (length <= L) up to just past supp g, coarse (length <= 1/2) beyond
        xm = min(self.xhi, self.supp[1] + Fraction(1, 10))
        nb1 = int((xm - self.xlo) / self.L) + 1
        grid = [self.xlo + k * (xm - self.xlo) / nb1 for k in range(nb1 + 1)]
        if xm < self.xhi:
            nb2 = int((self.xhi - xm) * 2) + 1
            grid += [xm + k * (self.xhi - xm) / nb2 for k in range(1, nb2 + 1)]
        nb = len(grid) - 1
        bps |= set(grid)
        bps = sorted(bps)
        nodeset = set(self.e)
        rootwin = [(lo, hi) for la, lo, hi, ball in self.R]
        import bisect
        efl = self.ef
        checks = 0; fails = 0; minlo = None; minat = None
        gi = 0
        for bk in range(nb):
            B0, B1 = grid[bk], grid[bk + 1]
            c = (B0 + B1) / 2; r = (B1 - B0) / 2
            cf, rf = fl(c), fl(r)
            near = [n for n in range(len(efl)) if abs(efl[n] - cf) <= 4 * rf + 1e-12 * 0] 
            # exact far test: |e_n - c| > 4r
            near = [n for n in range(bisect.bisect_left(efl, cf - 4.5 * rf), bisect.bisect_right(efl, cf + 4.5 * rf)) if abs(self.e[n] - c) <= 4 * r]
            nearset = set(near)
            far = [n for n in range(len(efl)) if n not in nearset]
            a, rem = self.block_taylor(c, r, far)
            ca = toarb(c)
            # cells in this block
            i0 = bisect.bisect_left(bps, B0); i1 = bisect.bisect_left(bps, B1)
            cells = [(bps[i], bps[i + 1]) for i in range(i0, i1)]
            stack = [(u, v, 0) for (u, v) in cells]
            while stack:
                u, v, dep = stack.pop()
                checks += 1
                lo = self.cell_lower(u, v, a, ca, rem, near)
                if not (lo > 0):
                    # certified counterexample? value of F at the midpoint (rigorous ball)
                    fm = self.point_value((u + v) / 2, a, ca, rem, near)
                    if fm < 0:
                        fails += 1
                        if fails <= verbose_fail: print(f"   COUNTEREXAMPLE F({fl((u+v)/2):.12f}) = {fm}  (certified negative)", flush=True)
                        continue
                    if dep < maxdepth:
                        mid = (u + v) / 2
                        stack.append((u, mid, dep + 1)); stack.append((mid, v, dep + 1))
                    else:
                        fails += 1
                        if fails <= verbose_fail: print(f"   FAIL [{fl(u):.12f},{fl(v):.12f}] lower {lo}", flush=True)
                else:
                    lf = float(lo.lower())
                    if minlo is None or lf < minlo: minlo = lf; minat = fl(u)
            if bk % 25 == 0:
                print(f"  block {bk}/{nb} x={fl(B0):.4f} checks={checks} fails={fails} min-lower={minlo} ({time.time()-t0:.0f}s)", flush=True)
        ok = okT and fails == 0
        print(f"main range [{self.xlo},{self.xhi}]: {fails==0}  ({checks} interval checks, min lower bound {minlo} at x={minat})  {time.time()-t0:.0f}s", flush=True)
        return ok

    def point_value(self, x, a, ca, rem, near):
        """rigorous ball for F(x) at a point x (Fraction) of the block, x not a node or root."""
        xa = toarb(x); val = xa + self.const
        for la, rlo, rhi, ball in self.R:
            if rlo <= x <= rhi: return arb(10**6)
            val += -la * abs(xa - ball).log()
        y = xa - ca; Pv = arb(0)
        for k in range(self.P, -1, -1): Pv = Pv * y + a[k]
        val += -2 * Pv
        val = val + arb(0).union(2 * rem).union(-2 * rem)
        for n in near:
            Cn = self.Ca[n]
            if Cn.is_zero(): continue
            val += -2 * Cn * G(xa - self.ea[n])
        return val

    def cell_lower(self, u, v, a, ca, rem, near):
        """rigorous lower bound of F on [u,v] (no node in (u,v))."""
        w = v - u; m = (u + v) / 2
        ma = toarb(m); ua = toarb(u); va = toarb(v); Ja = ua.union(va)
        hw2 = toarb(w) * toarb(w) / 4   # (w/2)^2
        # Rest value and derivative at m, Rest'' over J
        val = ma + self.const; der = arb(1); curv = arb(0)
        tan_v = arb(0); tan_d = arb(0)    # tangent terms (convex, singular at endpoints)
        chord_u = arb(0); chord_v = arb(0)
        const_extra = arb(0)
        for la, rlo, rhi, ball in self.R:
            if rhi < u or rlo > v:            # root strictly outside closed cell
                t = ma - ball
                val += -la * abs(t).log(); der += -la / t
                tJ = Ja - ball
                curv += la / (tJ * tJ)
            elif (rlo >= v or rhi <= u):     # touches an endpoint only (exact rational root at endpoint, or enclosure abutting)
                # -la log|x-r| convex on (u,v): tangent at m
                t = ma - ball
                tan_v += -la * abs(t).log(); tan_d += -la / t
            else:                             # root window: crude constant bound
                sup = max(abs(toarb(v) - toarb(rlo)).upper(), abs(toarb(u) - toarb(rhi)).upper())
                sup = arb(sup) + toarb(rhi - rlo)
                const_extra += -la * sup.log()
        # far Taylor polynomial in y = x - c
        y = ma - ca; yJ = Ja - ca
        Pv = arb(0); Pd = arb(0); Pdd = arb(0)
        Pk = self.P
        for k in range(Pk, -1, -1):
            Pv = Pv * y + a[k]
        for k in range(Pk, 0, -1):
            Pd = Pd * y + a[k] * k
        for k in range(Pk, 1, -1):
            Pdd = Pdd * yJ + a[k] * (k * (k - 1))
        val += -2 * Pv; der += -2 * Pd; curv += -2 * Pdd
        # near nodes
        for n in near:
            Cn = self.Ca[n]
            if Cn.is_zero(): continue
            en = self.e[n]
            if en == u or en == v:
                t = ma - self.ea[n]
                if self.C[n] > 0:   # convex term -2C G : tangent at m
                    tan_v += -2 * Cn * G(t); tan_d += -2 * Cn * G1(t)
                else:               # concave: chord
                    chord_u += -2 * Cn * G(ua - self.ea[n]); chord_v += -2 * Cn * G(va - self.ea[n])
            else:
                t = ma - self.ea[n]
                val += -2 * Cn * G(t); der += -2 * Cn * G1(t)
                tJ = Ja - self.ea[n]
                curv += -2 * Cn * abs(tJ).log()
        kappa = curv.lower()                      # inf of Rest'' over J (rigorous lower end)
        corr = arb(0)
        if kappa < 0: corr = -arb(kappa) * hw2 / 2   # positive amount subtracted
        hw = toarb(w) / 2
        base = val + tan_v - corr - 2 * rem + const_extra
        lo_u = base - (der + tan_d) * hw + chord_u
        lo_v = base + (der + tan_d) * hw + chord_v
        return arb(min(lo_u.lower(), lo_v.lower()))

if __name__ == '__main__':
    path = sys.argv[1]
    ch = Checker(path)
    if len(sys.argv) > 2 and sys.argv[2] == 'lambda':
        ch.lam = Fraction(sys.argv[3]); print("lambda overridden:", ch.lam)
    ch.setup_energy()
    ok = ch.run()
    print("VERDICT:", "VERIFIED" if ok else "NOT VERIFIED", f"C_86 >= {ch.lam} = {fl(ch.lam):.12f}")