#!/usr/bin/env python3
r"""lbl3disp.py -- shared symbolic-dispersion library for rows 28/29/33
(LBL3SE / LBL3KP / LBL3VP).  mpmath ONLY.  No AMFlow, no Kira, no file reads.

===========================  THE SYMBOLIC FORM  ==============================
The equal-mass kite self-energy spectral density rho(w) = -Im Sigma_kite(w+i0)
(m^2=1, w = l^2/m^2) is derived in closed form from the discontinuity of the
exact-rational kite differential system (row 8 of d/dw M = A(w,d) M, A linear
in d; exact rational, validated max-diff 0 against the IBP samples; entries
vendored below as explicit rationals).

All epsilon^-2 master discontinuities vanish; the eps^-1 discs are
Im B^(0) = pi (w-1)/w for the T*B and B^2 masters only; the kite disc is
finite.  A0[8,8] = -1/w exactly, so phi = Im M_kite^(0) obeys the scalar ODE
    d phi/dw = -phi/w + F(w),
and threshold regularity phi(1+) = 0 kills the integration constant:
    rho(w) = -(1/w) * int_1^w  v F(v) dv          [NO boundary input at all]

The inhomogeneity F collapses analytically (this work):
  * bubble masters j in {1,7} (T*B and B^2):   F_17(v) = 2 pi ln(v-1)/v^2
  * (m,0,0)-sunrise pair j in {5,6}:           F_56(v) = -pi ln v/(v(v-1))
    [Im S_m00^(0)(v) = -pi((v^2-1)/(2v) - ln v), Im S_m00dot^(0) = pi(ln v - (v-1)/v),
     both verified to 25 digits against recorded AMFlow master Laurent data at w=5]
  * (m,m,m)-sunrise pair j in {2,3} (v>9 only): F_23(v) below.
The elementary part integrates in closed form (dilogarithm level):
    rho_elem(w) = -(pi/w) [ 2 ln(w-1) ln w + 3 Li2(1-w) ]        (exact, all w>1)
which reproduces the threshold model rho ~ -2 pi (w-1) ln(w-1) analytically.

The elliptic (Gamma_1(6) sunrise-cut) part, for v > 9, with
xm=(sqrt v - 1)^2, xp=(sqrt v + 1)^2:
    ImS(v)    = -(pi/v) int_4^xm sqrt((xm-s)(xp-s)) sqrt((s-4)/s) ds
    ImSdot(v) = +(pi/v) int_4^xm (v-1+s) sqrt((s-4)/s) / sqrt((xm-s)(xp-s)) ds
                [d/d(m3^2) of the cut phase space; the bracket collapses to
                -2(v-1+s) exactly -- cancellation-free single positive-integrand
                quadrature, replaces the (2 ImS + pi dI3/dv)/3 DE-block form
                which loses ~v/ln v digits]
    F_23(v)   = [2 ImS - (v+3) ImSdot] / (v (v-1)^2)
I3 is the 3-body equal-mass phase-space integral = the sunrise spectral
density (the "known Gamma_1(6) object"); its normalization -pi was verified
exactly against recorded AMFlow Laurent data.

Full density:  rho(w) = rho_elem(w) - (1/w) int_9^w v F_23(v) dv   (w > 9).

Dispersive integrals (I(s,t,m2) = (1/pi) int_1^inf rho(w) K(w) dw) are
computed SPLIT (each piece absolutely convergent; the large-w cancellation
between rho_elem and the elliptic part then never occurs numerically):
    I = (1/pi) int_1^inf rho_elem(w) K(w) dw
      - (1/pi) int_9^inf v F_23(v) G(v) dv,     G(v) = int_v^inf K(w)/w dw
(second line = Fubini swap of the iterated integral; both pieces absolutely
convergent since v F_23 G ~ ln^2 v/v^3).

Kernels:
  row 28:  K = Box1(s,t;m2,M^2=w), the one-loop massive box in the closed
           1-dimensional dilogarithmic form (two Feynman integrations done
           analytically, 4-root partial fractions; port of
           the closed dilogarithmic box with Vieta-stabilized
           small roots for huge w).
  row 29:  K_P(w) = [K(w)-K(1)]/(w-1)^2 - K'(1)/(w-1)  and
           K_1(w) = [K(w)-K(1)]/(w-1)  (QED-parent stub kernels).
=============================================================================
CHANGELOG:
  2026-07-03: y=ln w map for the
  large-w tail (GMoment._direct + integrate_disp), lifting the x=1/w log-endpoint
  lower bound on the rate (~0.34 d/node -> ~0.8 d/node); cancellation-free single-integrand ImSdot
  (imS_pair/F23). A1 = -2pi is now a THEOREM of rho_elem, not a constant fit.
"""
import time

import mpmath as mp


# ------------------------------------------------------------------ helpers
def agree_digits(a, b):
    d = mp.fabs(a - b)
    if d == 0:
        return mp.inf
    return -mp.log10(d / mp.fabs(b))


# ---------------------------------------------------- rho: elementary part
def rho_elem(w):
    """rho below the 3-particle cut and the elementary part above it:
    -(pi/w)[2 ln(w-1) ln w + 3 Li2(1-w)]  (exact closed form)."""
    w = mp.mpf(w)
    u = w - 1
    if u <= 0:
        return mp.mpf(0)      # threshold limit (rho ~ -2 pi u ln u -> 0)
    return -(mp.pi / w) * (2 * mp.log(u) * mp.log(w) + 3 * mp.polylog(2, -u))


# ------------------------------------------- elliptic source F_23 (v > 9)
def _gl_nodes(n):
    key = ('gl', n, mp.mp.prec)
    cache = _gl_nodes.cache
    if key not in cache:
        cache[key] = mp.gauss_quadrature(n, 'legendre')
    return cache[key]


_gl_nodes.cache = {}


def _quad_gl(f, a, b, n):
    nodes, wts = _gl_nodes(n)
    c = (b - a) / 2
    d = (b + a) / 2
    return c * mp.fsum(wt * f(c * x + d) for x, wt in zip(nodes, wts))


def _quad_gl2(f, a, b, n):
    """Same but f returns a pair; both integrated in one pass."""
    nodes, wts = _gl_nodes(n)
    c = (b - a) / 2
    d = (b + a) / 2
    s1 = mp.mpf(0)
    s2 = mp.mpf(0)
    for x, wt in zip(nodes, wts):
        r1, r2 = f(c * x + d)
        s1 += wt * r1
        s2 += wt * r2
    return c * s1, c * s2


def imS_pair(v, dps=None):
    """(ImS, ImSdot)(v>9): equal-mass sunrise 3-body-cut densities, each as a
    SINGLE positive-integrand quadrature (no cancellation at any v):
      ImS    = -(pi/v) int_4^xm sqrt((xm-s)(xp-s)) sqrt((s-4)/s) ds
      ImSdot = +(pi/v) int_4^xm (v-1+s) sqrt((s-4)/s) / sqrt((xm-s)(xp-s)) ds
    (ImSdot from d/d(m3^2) of the cut phase-space integral at m3^2=1; the
    bracket -(sqrt v -1)(xp-s)+(sqrt v +1)(xm-s) = -2(v-1+s) exactly.)
    Endpoint-regularized (y^2 = s-4, z^2 = xm-s), geometric panels: accurate
    uniformly in v.  Verified against the {2,3} DE-block relation
    ImSdot = (2 ImS + pi dI3/dv)/3 and recorded AMFlow Laurent data."""
    if dps is None:
        dps = mp.mp.dps
    old = mp.mp.dps
    mp.mp.dps = dps + 10
    try:
        v = mp.mpf(v)
        sw = mp.sqrt(v)
        xm = (sw - 1) ** 2
        xp = (sw + 1) ** 2
        n = int(0.9 * (dps + 10)) + 12       # endpoint-map regions (A, first C)
        # 2026-07-06: slope 0.4->0.9 cures the ~0.44*dps+17.7 d imS_pair/F23 plateau
        nb = max(12, int(0.9 * (dps + 10)) + 8)   # geometric ln-panels (B, C-chain)

        def f_pair(s):
            root = mp.sqrt((xm - s) * (xp - s))
            fac = mp.sqrt((s - 4) / s)
            return root * fac, (v - 1 + s) * fac / root

        I = mp.mpf(0)
        D = mp.mpf(0)
        smid = min(mp.mpf(16), (4 + xm) / 2)
        # region A: [4, smid], s = 4 + y^2  (kills the sqrt(s-4) endpoint)
        Ya = mp.sqrt(smid - 4)

        def fa(y):
            s = 4 + y * y
            r1, r2 = f_pair(s)
            return r1 * 2 * y, r2 * 2 * y

        i_, d_ = _quad_gl2(fa, mp.mpf(0), Ya, n)
        I += i_; D += d_
        # region C: [cmid, xm], s = xm - z^2 (kills 1/sqrt(xm-s)); geometric
        # z-panels (analyticity radius ~ 2 v^(1/4) from xp - s = 4 sqrt(v) + z^2)
        cmid = max(smid, xm / 2) if xm / 2 > smid else smid
        Zc = mp.sqrt(xm - cmid)

        def fc(z):
            s = xm - z * z
            fac = mp.sqrt((s - 4) / s)
            rootp = mp.sqrt(xp - s)
            return fac * rootp * 2 * z * z, (v - 1 + s) * fac * 2 / rootp

        if Zc > 0:
            z1 = min(2 * v ** mp.mpf('0.25'), Zc)
            i_, d_ = _quad_gl2(fc, mp.mpf(0), z1, n)
            I += i_; D += d_
            a_ = z1
            while a_ < Zc * (1 - mp.mpf(10) ** (-dps)):
                b_ = min(a_ * mp.e, Zc)
                i_, d_ = _quad_gl2(fc, a_, b_, nb)
                I += i_; D += d_
                a_ = b_
        # region B: [smid, cmid] geometric panels in s
        if cmid > smid:
            a = smid
            while a < cmid * (1 - mp.mpf(10) ** (-dps)):
                b = min(a * mp.e, cmid)
                i_, d_ = _quad_gl2(f_pair, a, b, nb)
                I += i_; D += d_
                a = b
        ImS = -mp.pi / v * I
        ImSdot = mp.pi / v * D
        res = (+ImS, +ImSdot)
    finally:
        mp.mp.dps = old
    return res


def F23(v, dps=None):
    """Elliptic part of the inhomogeneity (v > 9): the Gamma_1(6) sunrise-cut
    density fed through the exact DE row-8 entries."""
    v = mp.mpf(v)
    if v <= 9:
        return mp.mpf(0)
    ImS, ImSdot = imS_pair(v, dps)
    return (2 * ImS - (v + 3) * ImSdot) / (v * (v - 1) ** 2)


# --------------------------------------------------- cumulative Phi_23, rho
class Phi23:
    """Phi(v) = int_9^v x F23(x) dx by descending-sorted cumulative panels
    (geometric GL sub-panels; F23 evals cached)."""

    def __init__(self, dps):
        self.dps = dps
        self.cache = {}          # mpf v -> Phi(v)
        self.F = {}

    def _f(self, x):
        key = mp.nstr(x, self.dps + 8)
        r = self.F.get(key)
        if r is None:
            r = x * F23(x, self.dps)
            self.F[key] = r
        return r

    def _panel(self, a, b):
        """int_a^b x F23 dx, geometric sub-panels ratio <= e."""
        tot = mp.mpf(0)
        n = int(0.8 * self.dps) + 12
        x0 = a
        while x0 < b * (1 - mp.mpf(10) ** (-self.dps - 5)):
            x1 = min(x0 * mp.e, b)
            tot += _quad_gl(self._f, x0, x1, n)
            x0 = x1
        return tot

    def value(self, v):
        v = mp.mpf(v)
        if v <= 9:
            return mp.mpf(0)
        if v in self.cache:
            return self.cache[v]
        # nearest cached value below/above in log distance
        best = None
        for u in self.cache:
            if best is None or abs(mp.log(v / u)) < abs(mp.log(v / best)):
                best = u
        if best is None:
            val = self._panel(mp.mpf(9), v)
        elif best < v:
            val = self.cache[best] + self._panel(best, v)
        else:
            val = self.cache[best] - self._panel(v, best)
        self.cache[v] = val
        return val


def make_rho(dps):
    """Full spectral density rho(w) (closed form + cumulative elliptic part)."""
    phi = Phi23(dps)

    def rho(w):
        w = mp.mpf(w)
        if w <= 1:
            return mp.mpf(0)
        r = rho_elem(w)
        if w > 9:
            r -= phi.value(w) / w
        return r
    rho.phi = phi
    return rho


# ------------------------------------------------------- Box1 kernel (dilog)
def box1(s, t, m2, M5, dps=40):
    """eps^0 of the one-loop massive box (pySecDec normalization), closed
    1D dilogarithmic form; deep-Euclidean s<0, t<0, u=-s-t<4 m2, M5>0.
    Vieta-stabilized small roots (mal = r1/(M5*malp)) so the rule is
    accurate for arbitrarily large M5."""
    work = dps + 18
    old = mp.mp.dps
    mp.mp.dps = work
    try:
        s = mp.mpf(s); t = mp.mpf(t); m2 = mp.mpf(m2); M5 = mp.mpf(M5)
        u = -s - t
        a = M5 + m2 - s
        if not (s < 0 and t < 0 and m2 > 0 and M5 > 0) or u >= 4 * m2:
            raise ValueError("box1: need s<0, t<0, m2>0, M5>0, u=-s-t<4m2")
        n = work + 12
        n += n % 2   # even rule: an odd rule places a node at the midpoint v = 2 exactly, where
                     # (u + s*malp) and (u + s*mbep) vanish together as M5 -> inf at kinematic
                     # points with -u/s = 3 (t = 2s), a removable singularity that the node would
                     # divide by; an even rule has no node there.
        nodes, weights = _gl_nodes(n)
        L = mp.mpf(2)
        one = mp.mpf(1)
        tot = mp.mpf(0)
        for tk, wk in zip(nodes, weights):
            vv = L * (one + tk) / (one - tk)
            dv = L * 2 / (one - tk) ** 2
            vp1 = vv + 1
            p1 = a * vv + (M5 + m2)
            r1 = m2 * vp1 * vp1
            p2 = (M5 + m2) * vp1
            r2 = r1 - u * vv
            sd1 = mp.sqrt(p1 * p1 - 4 * M5 * r1)
            sd2 = mp.sqrt(p2 * p2 - 4 * M5 * r2)
            twoM5 = 2 * M5
            malp = (p1 + sd1) / twoM5
            mal = r1 / (M5 * malp) if malp != 0 else (p1 - sd1) / twoM5
            mbep = (p2 + sd2) / twoM5
            mbe = r2 / (M5 * mbep) if mbep != 0 else (p2 - sd2) / twoM5
            gv = -((-1 / ((u + s * mal) * sd1)) * mp.log(mal)
                   + (1 / ((u + s * malp) * sd1)) * mp.log(malp)
                   + (1 / ((u + s * mbe) * sd2)) * mp.log(mbe)
                   + (-1 / ((u + s * mbep) * sd2)) * mp.log(mbep))
            tot += wk * gv * dv
        res = +tot
    finally:
        mp.mp.dps = old
    return res


def make_K(s, t, m2, dps):
    cache = {}

    def K(w):
        key = mp.nstr(w, dps + 8)
        r = cache.get(key)
        if r is None:
            r = box1(s, t, m2, w, dps)
            cache[key] = r
        return r
    K.cache = cache
    return K


# --------------------------------------- kernel log-moment G(v)=int K/w dw
class GMoment:
    """G(v) = int_v^inf kern(w)/w dw, cumulative with a direct mapped-TS
    reference value at the largest v ever seen."""

    def __init__(self, kern, dps):
        self.kern = kern
        self.dps = dps
        self.cache = {}

    def _direct(self, v):
        """int_v^inf kern(w)/w dw = int_{ln v}^inf kern(e^y) dy: in y the
        integrand is analytic (log structure -> polynomial) and ~ e^{-2y} poly;
        geometric GL y-panels, stop on the remaining-integral bound."""
        import math
        thresh = mp.mpf(10) ** (-(self.dps + 9))
        Ly = mp.mpf(6)
        n = int((self.dps + 9) / 0.8) + 8

        def f(y):
            return self.kern(mp.exp(y))

        y0 = mp.log(v)
        tot = mp.mpf(0)
        npan = 0
        while True:
            tot += _quad_gl(f, y0, y0 + Ly, n)
            y0 += Ly
            npan += 1
            if 1000 * (y0 + 1) ** 3 * mp.exp(-y0) < thresh or npan > 60 * self.dps:
                break
        return tot

    def _panel(self, a, b):
        """int_a^b kern/w dw = int kern d(ln w), geometric sub-panels."""
        tot = mp.mpf(0)
        x0 = a
        while x0 < b * (1 - mp.mpf(10) ** (-self.dps - 5)):
            x1 = min(x0 * mp.e, b)
            Llen = float(mp.log(x1 / x0))
            # GL error ~ (Llen/(2 rho))^(2n); guarantee < 1e-(dps+8) for any Llen
            import math
            n = max(8, int(2.303 * (self.dps + 8) / (2 * max(1.5, math.log(8.0 / max(Llen, 1e-12))))) + 4)
            tot += _quad_gl(lambda w: self.kern(w) / w, x0, x1, n)
            x0 = x1
        return tot

    def value(self, v):
        v = mp.mpf(v)
        if v in self.cache:
            return self.cache[v]
        best = None
        for u in self.cache:
            if best is None or abs(mp.log(v / u)) < abs(mp.log(v / best)):
                best = u
        if best is None:
            val = self._direct(v)
        elif best > v:
            val = self.cache[best] + self._panel(v, best)
        else:
            val = self.cache[best] - self._panel(best, v)
        self.cache[v] = val
        return val


# ------------------------------------------------------ TS panel machinery
def ts_ladder(f, a, b, dps, max_level=9, tol_digits=None):
    """Tanh-sinh ladder on [a,b]; returns (best, ladder list). Stops when two
    consecutive levels agree to tol."""
    old = mp.mp.dps
    mp.mp.dps = dps + 10
    try:
        ts = mp.calculus.quadrature.TanhSinh(mp.mp)
        prec = mp.mp.prec
        tol = mp.mpf(10) ** (-(tol_digits or (dps + 6)))
        results = []
        for m in range(1, max_level + 1):
            nodes = ts.get_nodes(mp.mpf(a), mp.mpf(b), m, prec)
            results.append(ts.sum_next(f, nodes, m, prec, results, False))
            if len(results) >= 3:
                d1 = abs(results[-1] - results[-2])
                if d1 <= tol * max(mp.mpf(1), abs(results[-1])):
                    break
        res = +results[-1]
        lad = [+x for x in results]
    finally:
        mp.mp.dps = old
    return res, lad


def integrate_disp(fw, dps, panels=(1, 2, 9, 200), tail_from=200, max_level=8,
                   methods=None, tail_bound=None):
    """int fw(w) dw over [panels...] + mapped [tail_from, inf).
    methods: per-panel 'ts' or 'gl' (default: 'ts' for first panel, 'gl' rest).
    'gl' panels are split geometrically (ln-ratio <= 1.2) with n ~ 1.1 dps.
    tail_bound(w): cheap upper bound on |fw(w)|; tail nodes with bound below
    1e-(dps+8) * scale are skipped (avoids astronomically deep evaluations)."""
    tot = mp.mpf(0)
    diags = []
    if methods is None:
        methods = ['ts'] + ['gl'] * (len(panels) - 2)
    ngl = int(1.1 * dps) + 12
    for (a, b), meth in zip(zip(panels[:-1], panels[1:]), methods):
        a = mp.mpf(a); b = mp.mpf(b)
        if meth == 'ts':
            r, lad = ts_ladder(fw, a, b, dps, max_level)
            tot += r
            diags.append((float(a), float(b), 'ts', len(lad),
                          float(agree_digits(lad[-1], lad[-2])) if len(lad) > 1 else None))
        else:
            r = mp.mpf(0)
            x0 = a
            while x0 < b * (1 - mp.mpf(10) ** (-dps - 5)):
                x1 = min(x0 * mp.exp(mp.mpf('1.2')), b)
                r += _quad_gl(fw, x0, x1, ngl)
                x0 = x1
            tot += r
            diags.append((float(a), float(b), 'gl', ngl, None))
    # tail: w = e^y, int_{tail_from}^inf f(w) dw = int_{y0}^inf f(e^y) e^y dy.
    # In y the integrand is ANALYTIC in a strip (all ln w -> polynomial in y)
    # and decays like e^{-y} f(e^y); GL panels of fixed y-length resolve it at
    # ~0.8 digits/node.  Stop when the panel bound falls below threshold.
    Wc = mp.mpf(tail_from)
    thresh = mp.mpf(10) ** (-(dps + 9))
    Ctail = mp.mpf(tail_bound if tail_bound is not None else 1000)

    def ftail(y):
        w = mp.exp(y)
        return fw(w) * w

    Ly = mp.mpf(6)
    ntail = int((dps + 9) / 0.8) + 8
    y0 = mp.log(Wc)
    rtail = mp.mpf(0)
    npan = 0
    while True:
        y1 = y0 + Ly
        rtail += _quad_gl(ftail, y0, y1, ntail)
        y0 = y1
        npan += 1
        # remaining-integral bound: |f(w)| <= Ctail ln^3(w)/w^2 assumed =>
        # int_{e^y0}^inf <= Ctail (y0+1)^3 e^{-y0}
        if Ctail * (y0 + 1) ** 3 * mp.exp(-y0) < thresh or npan > 60 * dps:
            break
    tot += rtail
    diags.append(('tail', 'inf', 'gl-y', npan, ntail))
    return tot, diags


# --------------------------------------------------------- full assemblies
def assemble(kern_K, Gm, dps, verbose=True, tail_from=200):
    """I = (1/pi) int rho_elem K dw - (1/pi) int_9^inf v F23(v) G(v) dv."""
    t0 = time.time()

    Ia, dga = integrate_disp(lambda w: rho_elem(w) * kern_K(w), dps,
                             panels=(1, 2, 9, tail_from), tail_from=tail_from,
                             tail_bound=1000)
    ta = time.time() - t0
    if verbose:
        print(f"    [Ia elementary x kernel] {mp.nstr(Ia / mp.pi, 25)}  ({ta:.1f}s; panels {dga})")
    t0 = time.time()
    Fc = {}

    def fb(v):
        key = mp.nstr(v, dps + 8)
        r = Fc.get(key)
        if r is None:
            r = v * F23(v, dps)
            Fc[key] = r
        return r * Gm.value(v)

    Ib, dgb = integrate_disp(fb, dps, panels=(9, 30, tail_from), tail_from=tail_from,
                             methods=['gl', 'gl'], tail_bound=1000)
    tb = time.time() - t0
    if verbose:
        print(f"    [Ib elliptic x G-moment] {mp.nstr(-Ib / mp.pi, 25)}  ({tb:.1f}s; panels {dgb})")
    return (Ia - Ib) / mp.pi


def assemble_control(dps, verbose=True, tail_from=200):
    """Kernel-swap control: (1/pi) int rho(w)/(w+1) dw  ==  -Sigma_kite(-1).
    G-tilde(v) = ln(1+1/v) EXACT."""
    Ia, _ = integrate_disp(lambda w: rho_elem(w) / (w + 1), dps,
                           panels=(1, 2, 9, tail_from), tail_from=tail_from,
                           tail_bound=1000)
    Fc = {}

    def fb(v):
        key = mp.nstr(v, dps + 8)
        r = Fc.get(key)
        if r is None:
            r = v * F23(v, dps)
            Fc[key] = r
        return r * mp.log1p(1 / v)

    Ib, _ = integrate_disp(fb, dps, panels=(9, 30, tail_from), tail_from=tail_from,
                           methods=['gl', 'gl'], tail_bound=1000)
    return (Ia - Ib) / mp.pi