#!/usr/bin/env python3
r"""k33lib.py (vendored subset) -- Gauss-Legendre node cache and the cancellation-free power
difference used by the row-33 eps^0 engine (row33_eps0.py).

Two functions of the research module k33lib.py, byte for byte (the rest of that module -- the
fixed-eps Feynman-parameter masters of the one-loop box family -- is not needed by this leg):
  gl_nodes(n, work)        Gauss-Legendre nodes/weights at `work` digits, cached per (n, work)
  _powdiff(A, B, X, p)     [A^{-p} - (A+B*X)^{-p}] / B, stable as B*X/A -> 0 (log1p / expm1 branch)
"""
import mpmath as mp

_GL_CACHE = {}


def gl_nodes(n, work):
    key = (n, work)
    xw = _GL_CACHE.get(key)
    if xw is None:
        old = mp.mp.dps
        mp.mp.dps = work
        xw = mp.gauss_quadrature(n, 'legendre')
        mp.mp.dps = old
        _GL_CACHE[key] = xw
    return xw


def _powdiff(A, B, X, p):
    """[A^{-p} - (A+B*X)^{-p}] / B  , stable as B*X/A -> 0 (p>0)."""
    z = B * X / A
    if abs(z) < mp.mpf('0.5'):
        # A^{-p} * (1 - (1+z)^{-p}) / B  = A^{-p} * (-expm1(-p*log1p(z))) / B
        if B == 0:
            return p * X * A ** (-p - 1)
        return A ** (-p) * (-mp.expm1(-p * mp.log1p(z))) / B
    return (A ** (-p) - (A + B * X) ** (-p)) / B