#!/usr/bin/env python3
r"""lbl3disp.py (vendored subset) -- the eps^0 one-loop massive box kernel in its closed
dilogarithmic form, used by the row-33 eps^0 engine (row33_eps0.py) as an INDEPENDENT check of the
transported kernel layer K^(0)(w).

One function of the shared research module lbl3disp.py, byte for byte, with its node-cache helper
(the rest of that module -- the kite spectral density and the dispersive assemblies of the sibling
rows -- is not needed by this leg):
  box1(s, t, m2, M5, dps)  eps^0 of the one-loop box with three equal masses m2 and one mass M5,
                           deep-Euclidean s<0, t<0, u=-s-t<4 m2; two Feynman integrations done
                           analytically, one Gauss-Legendre fold; Vieta-stabilized small roots so the
                           rule stays accurate for arbitrarily large M5.
mpmath only.
"""
import mpmath as mp


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 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