#!/usr/bin/env python3 r"""Row 33 kernel de-AMFlow: exact fixed-eps Feynman-parameter evaluators for the one-loop box family (prop masses (wp,1,1,1), massless legs; kinematics (s,u) are PARAMETERS since build B1, 2026-09-04, defaults = the reference s=-1, u=4/3). Family: D1 = l^2 - wp, D2 = (l+p1)^2 - 1, D3 = (l+p1+p2)^2 - 1, D4 = (l-p4)^2 - 1 p_i^2 = 0, (p1+p2)^2 = s, (p2+p3)^2 = u [t enters via u=-s-t; reference s=-1, u=4/3, t=-1/3; every function below takes (s, u) as EXACT rationals (int / fractions.Fraction / sympy.Rational), converted to mpf inside the function at the caller's precision -- at the defaults the arithmetic is the record's, bit for bit: -s = mpf(1) multiplies exactly, u = mpf(4)/3 as before] Normalization (VERIFIED vs the recorded 1-loop AMFlow seeds: plain measure, no e^{gamma eps}): I_N = (-1)^N Gamma(N-d/2) int_simplex F^{d/2-N}, d = 4-2eps F = sum x_i m_i^2 - sum_{i0 everywhere here: Euclidean) Masters and their F on the unit simplex (x's ordered by the props PRESENT); general (s,u), with q_13^2 = s (D1-D3 carry p1+p2) and q_24^2 = u (D2-D4 carry p1+p4 = -(p2+p3)), all other q_ij^2 = 0: T(m2) = Gamma(eps)/(1-eps) * m2^{1-eps} [closed] bub_s (1,0,1,0): F = 1 + x(wp-1) - s*x(1-x) (q2=s) [1-fold] bub_u (0,1,0,1): F = 1 - u*x(1-x) (q2=u) [1-fold, wp-free] tri_s (1,1,1,0): F = 1 + x1(wp-1) - s*x1*x3 (x3 folded) [1-fold] tri_u (1,1,0,1): F = 1 + x1(wp-1) - u*x2*x4 (x1 folded) [1-fold] tri_c (0,1,1,1): F = 1 - u*xa*xc (xa folded) [1-fold, wp-free] box (1,1,1,1): F = 1 + x1(wp-1) - s*x1*x3 - u*x2*x4 x3-fold analytic: F = A + B*x3, B = -s*x1 + u*x2, A = 1 + x1(wp-1) - u*x2*(1-x1-x2), A+B*X = 1+x1(wp-1)-s*x1*X, X = 1-x1-x2 [2-fold] (at s=-1, u=4/3 these are the record's forms: +x1*x3, -(4/3)x2*x4, B=x1+(4/3)x2) Kn threshold Taylor (wp = 1+u): F = F1 + u*x1 => K_n = Gamma(2+eps) * (2+eps)_n/n! * (-1)^n * int x1^n F1^{-2-eps-n} (EXACT moments, same x3-fold, shared node grid) All integrands analytic on the closed domains for wp > ~0.9 (F >= 1-u/4 - ish), so plain Gauss-Legendre converges geometrically (rate calibrated at u=4/3; degrades as u -> 4). Stable small-z branch via log1p/expm1. CHANGELOG: 2026-09-04 build B1 (the (s,t)-general build): (s, u) parameters on bub_s, bub_u, tri_s, tri_u, tri_c, box, kn_moments, masters_at; defaults S_REF=-1, U_REF=4/3 reproduce the 2026-07-03 file bit for bit (regression-checked). """ from fractions import Fraction import mpmath as mp _GL_CACHE = {} S_REF = Fraction(-1) # reference (p1+p2)^2 U_REF = Fraction(4, 3) # reference (p2+p3)^2 = -s-t at t=-1/3 def _mpq(q): """Exact rational (int / Fraction / sympy.Rational / 'p/q' string) -> mpf at the CURRENT precision, as mpf(p)/q (the record's `mp.mpf(4) / 3` form).""" q = Fraction(q) if not isinstance(q, Fraction) else q return mp.mpf(q.numerator) / q.denominator def _ms(s): """-s as mpf (the coefficient of x1*x3 in F); mpf(1) at the reference.""" return _mpq(-Fraction(s)) 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 # ---------------- closed forms ---------------- def tadpole(m2, eps): return mp.gamma(eps) / (1 - eps) * mp.mpf(m2) ** (1 - eps) # ---------------- 1-fold masters ---------------- def _quad01(f, n, work): nodes, weights = gl_nodes(n, work) half = mp.mpf('0.5') tot = mp.mpf(0) for t, w in zip(nodes, weights): x = half * (t + 1) tot += w * f(x) return half * tot def bub_s(wp, eps, n, work, s=S_REF): """(1,0,1,0): masses (wp,1) at q^2=s. Gamma(eps) int_0^1 F^-eps, F = 1 + x(wp-1) - s*x(1-x).""" wp = mp.mpf(wp) ms = _ms(s) f = lambda x: (1 + x * (wp - 1) + ms * x * (1 - x)) ** (-eps) return mp.gamma(eps) * _quad01(f, n, work) def bub_u(eps, n, work, u=U_REF): """(0,1,0,1): masses (1,1) at q^2=u. F = 1 - u*x(1-x).""" c = _mpq(u) f = lambda x: (1 - c * x * (1 - x)) ** (-eps) return mp.gamma(eps) * _quad01(f, n, work) def tri_s(wp, eps, n, work, s=S_REF): """(1,1,1,0): F = 1 + x1(wp-1) - s*x1*x3, x3 in [0,1-x1] folded: int dx3 (A+C x3)^{-1-eps} = powdiff(A,C,X,eps)/eps with C=-s*x1.""" wp = mp.mpf(wp) ms = _ms(s) def f(x1): A = 1 + x1 * (wp - 1) return _powdiff(A, ms * x1, 1 - x1, eps) / eps return -mp.gamma(1 + eps) * _quad01(f, n, work) def tri_u(wp, eps, n, work, u=U_REF): """(1,1,0,1): F = 1 + x1(wp-1) - u*x2*x4, x4=1-x1-x2. Fold x1 over [0,1-x2]: F = A + C*x1, A = 1-u*x2(1-x2), C = wp-1+u*x2.""" wp = mp.mpf(wp) c = _mpq(u) def f(x2): A = 1 - c * x2 * (1 - x2) C = wp - 1 + c * x2 return _powdiff(A, C, 1 - x2, eps) / eps return -mp.gamma(1 + eps) * _quad01(f, n, work) def tri_c(eps, n, work, u=U_REF): """(0,1,1,1): F = 1 - u*xa*xc. Fold xa over [0,1-xc]: F = 1 + C*xa, C = -u*xc.""" c = _mpq(u) def f(xc): return _powdiff(mp.mpf(1), -c * xc, 1 - xc, eps) / eps return -mp.gamma(1 + eps) * _quad01(f, n, work) # ---------------- box (2-fold) ---------------- def box(wp, eps, n, work, s=S_REF, u=U_REF): """(1,1,1,1): after the analytic x3-fold, box = Gamma(2+eps)/(1+eps) * int_tri dx1 dx2 powdiff(A, B, X, 1+eps) A = 1+x1(wp-1)-u*x2*X, B = -s*x1+u*x2, X = 1-x1-x2.""" wp = mp.mpf(wp) c = _mpq(u) ms = _ms(s) nodes, weights = gl_nodes(n, work) half = mp.mpf('0.5') p = 1 + eps tot = mp.mpf(0) for t1, w1 in zip(nodes, weights): x1 = half * (t1 + 1) L = 1 - x1 # x2 in [0, L] acc = mp.mpf(0) for t2, w2 in zip(nodes, weights): x2 = L * half * (t2 + 1) X = 1 - x1 - x2 A = 1 + x1 * (wp - 1) - c * x2 * X B = ms * x1 + c * x2 acc += w2 * _powdiff(A, B, X, p) tot += w1 * acc * (L * half) return mp.gamma(2 + eps) / (1 + eps) * tot * half def kn_moments(eps, nmax, n, work, s=S_REF, u=U_REF): """K_j(eps), j=0..nmax-1: exact Taylor coefficients of box(1+u) in u (here u is the threshold offset wp-1; the kinematic u is the parameter). box(1+u) = Gamma(2+eps) int (F1 + u*x1)^{-2-eps}, K_j = Gamma(2+eps) * (-1)^j * binom(-2-eps, j)... explicitly = Gamma(2+eps) * (2+eps)_j / j! * (-1)^j * M_j, M_j = int x1^j F1^{-2-eps-j} (sign: (F1+u x1)^{-2-eps} = sum_j C(-2-eps,j) u^j x1^j F1^{-2-eps-j}, C(-2-eps,j) = (-1)^j (2+eps)_j / j!) Each M_j gets its own analytic x3-fold: F1 = A1 + B*x3 => M_j = int_tri x1^j powdiff(A1, B, X, 1+eps+j)/(1+eps+j). All j share one (x1,x2) grid; powers built recursively per node.""" c = _mpq(u) ms = _ms(s) nodes, weights = gl_nodes(n, work) half = mp.mpf('0.5') tots = [mp.mpf(0)] * nmax for t1, w1 in zip(nodes, weights): x1 = half * (t1 + 1) L = 1 - x1 accs = [mp.mpf(0)] * nmax for t2, w2 in zip(nodes, weights): x2 = L * half * (t2 + 1) X = 1 - x1 - x2 A1 = 1 - c * x2 * X # F(wp=1, x3=0) B = ms * x1 + c * x2 x1p = mp.mpf(1) for j in range(nmax): accs[j] += w2 * x1p * _powdiff(A1, B, X, 1 + eps + j) x1p *= x1 for j in range(nmax): tots[j] += w1 * accs[j] * (L * half) G = mp.gamma(2 + eps) out = [] poch = mp.mpf(1) # (2+eps)_j / j! for j in range(nmax): if j > 0: poch *= (1 + eps + j) / j out.append(G * poch * (-1) ** j * tots[j] * half / (1 + eps + j)) return out def masters_at(wp, eps, n, work, s=S_REF, u=U_REF): """All 8 one-loop box masters at fixed eps and kinematics (s, u), in the reference master order: [1111, 1000, 0100, 1010, 0101, 1110, 1101, 0111].""" return [ box(wp, eps, n, work, s, u), tadpole(wp, eps), tadpole(1, eps), bub_s(wp, eps, n, work, s), bub_u(eps, n, work, u), tri_s(wp, eps, n, work, s), tri_u(wp, eps, n, work, u), tri_c(eps, n, work, u), ]