Notebooks · posts/the-pion-loop-function-by-computer-algebra

The two-pion loop function, built and differentiated by a small computer algebra system

Adam Zettel · September 25, 2026 · written with an agent

The loop function $L (q)$ sets the momentum dependence of the two-pion-exchange nuclear force at next-to-leading order. The program writes it as an expression tape, evaluates it, and differentiates the tape symbolically into a second tape. Both are checked against a 30-digit Feynman-parameter integral that never uses the closed form: $L (200 MeV) = 1.14642288565$ and $d L / d q = 1.24075106 \times 1 0^{- 3}$ /MeV there, to $1 0^{- 14}$ and $2 \times 1 0^{- 17}$ /MeV.

L (q) = \frac{w}{q} ln \frac{w + q}{2 m _{π}}, w = 4 m_{π}^{2} + q^{2}; L (q) = 1 + \frac{1}{2} \int_{0}^{1} ln (1 + x (1 - x) \frac{q ^{2}}{m _{π}^{2}}) d x

In chiral effective field theory the force between two nucleons is a sum of pion exchanges and short-range contact terms. At next-to-leading order the exchange of two pions enters through one function of the momentum transfer $q$ ,

L (q) = \frac{w}{q} ln \frac{w + q}{2 m _{π}}, w = 4 m_{π}^{2} + q^{2},

which multiplies the central, tensor and spin-spin parts of the potential (Machleidt and Entem, Phys. Rep. 503, 1 (2011), Eqs. (4.9)-(4.12), with their average pion mass $m_{π} = 138.039$ MeV). $L (0) = 1$ , and $L$ grows like $ln q$ at large $q$ . As a function of $q^{2}$ it has a branch cut that starts at $q^{2} = - 4 m_{π}^{2}$ , the two-pion threshold, so in coordinate space the force it carries has range $1/ (2 m_{π})$ , half the one-pion range.

Higher orders need many such functions and their derivatives, which is the job of a computer algebra system. This post builds a small one: an expression tape, a list of operations each of which refers to earlier entries, with an evaluator and a differentiator that turns a tape into the tape of its derivative.

The reference, posts/refs/the-pion-loop-function-by-computer-algebra.py, does not use the closed form. It computes the loop from its Feynman-parameter integral (the equation under the title) with mpmath's quadrature at 30 digits, and the derivative by mpmath's numerical differentiation of that integral.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 38–156

MPI = 138.039
# the most entries a tape can hold
TAPE_MAX = 512

# opcodes
OP_CONST = 0
OP_VAR = 1
OP_ADD = 2
OP_SUB = 3
OP_MUL = 4
OP_DIV = 5
OP_SQRT = 6
OP_LOG = 7

# the tape: entry k is op[k] applied to entries arg_a[k] and arg_b[k] (or holds val[k])
t_op = zeros(TAPE_MAX)
t_a = zeros(TAPE_MAX)
t_b = zeros(TAPE_MAX)
t_val = zeros(TAPE_MAX)
t_len = zeros(1)
# work space: the value of each entry, and the derivative entry of each entry
t_ev = zeros(TAPE_MAX)
t_dd = zeros(TAPE_MAX)

#| Append an entry to the tape; returns its index.
def push(o, a, b, c) {
    i = t_len[0]
    t_op[i] = o
    t_a[i] = a
    t_b[i] = b
    t_val[i] = c
    t_len[0] = i + 1
    return i
}

def cst(v) {
    return push(OP_CONST, 0, 0, v)
}
def var() {
    return push(OP_VAR, 0, 0, 0.0)
}
def add(a, b) {
    return push(OP_ADD, a, b, 0.0)
}
def sub(a, b) {
    return push(OP_SUB, a, b, 0.0)
}
def mul(a, b) {
    return push(OP_MUL, a, b, 0.0)
}
def div(a, b) {
    return push(OP_DIV, a, b, 0.0)
}
def tsqrt(a) {
    return push(OP_SQRT, a, 0, 0.0)
}
def tlog(a) {
    return push(OP_LOG, a, 0, 0.0)
}

#| Evaluate the tape at x in one forward pass; returns the value of entry `root`.
def eval_at(root, x) {
    for k to t_len[0] {
        o = t_op[k]
        a = t_a[k]
        b = t_b[k]
        if o == OP_CONST {
            t_ev[k] = t_val[k]
        } else if o == OP_VAR {
            t_ev[k] = x
        } else if o == OP_ADD {
            t_ev[k] = t_ev[a] + t_ev[b]
        } else if o == OP_SUB {
            t_ev[k] = t_ev[a] - t_ev[b]
        } else if o == OP_MUL {
            t_ev[k] = t_ev[a] * t_ev[b]
        } else if o == OP_DIV {
            t_ev[k] = t_ev[a] / t_ev[b]
        } else if o == OP_SQRT {
            t_ev[k] = sqrt(t_ev[a])
        } else {
            t_ev[k] = log(t_ev[a])
        }
    }
    return t_ev[root]
}

#| Symbolic differentiation, tape to tape: for each entry k in order, append
#| entries that compute its derivative and record where, in t_dd[k]. The
#| rules are the sum, product and quotient rules, and
#| sqrt(u)' = u' / (2 sqrt(u)) and ln(u)' = u' / u; the derivative refers to
#| the original entries, so both live on one tape. Returns the entry of
#| d(root)/dx.
def diff(root) {
    n = t_len[0]
    for k to n {
        o = t_op[k]
        a = t_a[k]
        b = t_b[k]
        if o == OP_CONST {
            t_dd[k] = cst(0.0)
        } else if o == OP_VAR {
            t_dd[k] = cst(1.0)
        } else if o == OP_ADD {
            t_dd[k] = add(t_dd[a], t_dd[b])
        } else if o == OP_SUB {
            t_dd[k] = sub(t_dd[a], t_dd[b])
        } else if o == OP_MUL {
            t_dd[k] = add(mul(t_dd[a], b), mul(a, t_dd[b]))
        } else if o == OP_DIV {
            t_dd[k] = div(sub(mul(t_dd[a], b), mul(a, t_dd[b])), mul(b, b))
        } else if o == OP_SQRT {
            t_dd[k] = div(t_dd[a], mul(cst(2.0), k))
        } else {
            t_dd[k] = div(t_dd[a], a)
        }
    }
    return t_dd[root]
}

MPI = 138.039
# the most entries a tape can hold
TAPE_MAX = 512

# opcodes
OP_CONST = 0
OP_VAR = 1
OP_ADD = 2
OP_SUB = 3
OP_MUL = 4
OP_DIV = 5
OP_SQRT = 6
OP_LOG = 7

# the tape: entry k is op[k] applied to entries arg_a[k] and arg_b[k] (or holds val[k])
t_op = zeros(TAPE_MAX)
t_a = zeros(TAPE_MAX)
t_b = zeros(TAPE_MAX)
t_val = zeros(TAPE_MAX)
t_len = zeros(1)
# work space: the value of each entry, and the derivative entry of each entry
t_ev = zeros(TAPE_MAX)
t_dd = zeros(TAPE_MAX)

#| Append an entry to the tape; returns its index.
def push(o, a, b, c) {
    i = t_len[0]
    t_op[i] = o
    t_a[i] = a
    t_b[i] = b
    t_val[i] = c
    t_len[0] = i + 1
    return i
}

def cst(v) {
    return push(OP_CONST, 0, 0, v)
}
def var() {
    return push(OP_VAR, 0, 0, 0.0)
}
def add(a, b) {
    return push(OP_ADD, a, b, 0.0)
}
def sub(a, b) {
    return push(OP_SUB, a, b, 0.0)
}
def mul(a, b) {
    return push(OP_MUL, a, b, 0.0)
}
def div(a, b) {
    return push(OP_DIV, a, b, 0.0)
}
def tsqrt(a) {
    return push(OP_SQRT, a, 0, 0.0)
}
def tlog(a) {
    return push(OP_LOG, a, 0, 0.0)
}

#| Evaluate the tape at x in one forward pass; returns the value of entry `root`.
def eval_at(root, x) {
    for k to t_len[0] {
        o = t_op[k]
        a = t_a[k]
        b = t_b[k]
        if o == OP_CONST {
            t_ev[k] = t_val[k]
        } else if o == OP_VAR {
            t_ev[k] = x
        } else if o == OP_ADD {
            t_ev[k] = t_ev[a] + t_ev[b]
        } else if o == OP_SUB {
            t_ev[k] = t_ev[a] - t_ev[b]
        } else if o == OP_MUL {
            t_ev[k] = t_ev[a] * t_ev[b]
        } else if o == OP_DIV {
            t_ev[k] = t_ev[a] / t_ev[b]
        } else if o == OP_SQRT {
            t_ev[k] = sqrt(t_ev[a])
        } else {
            t_ev[k] = log(t_ev[a])
        }
    }
    return t_ev[root]
}

#| Symbolic differentiation, tape to tape: for each entry k in order, append
#| entries that compute its derivative and record where, in t_dd[k]. The
#| rules are the sum, product and quotient rules, and
#| sqrt(u)' = u' / (2 sqrt(u)) and ln(u)' = u' / u; the derivative refers to
#| the original entries, so both live on one tape. Returns the entry of
#| d(root)/dx.
def diff(root) {
    n = t_len[0]
    for k to n {
        o = t_op[k]
        a = t_a[k]
        b = t_b[k]
        if o == OP_CONST {
            t_dd[k] = cst(0.0)
        } else if o == OP_VAR {
            t_dd[k] = cst(1.0)
        } else if o == OP_ADD {
            t_dd[k] = add(t_dd[a], t_dd[b])
        } else if o == OP_SUB {
            t_dd[k] = sub(t_dd[a], t_dd[b])
        } else if o == OP_MUL {
            t_dd[k] = add(mul(t_dd[a], b), mul(a, t_dd[b]))
        } else if o == OP_DIV {
            t_dd[k] = div(sub(mul(t_dd[a], b), mul(a, t_dd[b])), mul(b, b))
        } else if o == OP_SQRT {
            t_dd[k] = div(t_dd[a], mul(cst(2.0), k))
        } else {
            t_dd[k] = div(t_dd[a], a)
        }
    }
    return t_dd[root]
}

The tape of $L (q)$

The program writes $L$ the way it is printed: $w$ as the square root of $4 m_{π}^{2} + q^{2}$ , then $(w / q) ln ((w + q) / (2 m_{π}))$ . That is 11 entries. Differentiating appends the derivative of every entry, including the derivatives of constants, which are zero; the tape grows to 36 entries. Nothing simplifies the result, so the derivative tape does more arithmetic than a hand-written derivative would, but it is exact.

[2]lines 169–181

q = var()
w = tsqrt(add(cst(4.0 * MPI * MPI), mul(q, q)))
L_root = mul(div(w, q), tlog(div(add(w, q), cst(2.0 * MPI))))
n_L = t_len[0]
dL_root = diff(L_root)
n_all = t_len[0]
show n_L
show n_all

test "the tape of L has 11 entries, and 36 once differentiated" {
    assert_near(n_L, 11, 0.0)
    assert_near(n_all, 36, 0.0)
}

Values against the Feynman-parameter integral

The reference evaluates $L$ at eight momenta from its integral representation, which comes from the loop diagram itself, not from the closed form. The tape agrees at every one to $1 0^{- 14}$ , from $q = 1$ MeV, where $L - 1 = 4.4 \times 1 0^{- 6}$ , to 1 GeV. At small $q$ the logarithm in the integral expands, and with $t = q^{2} / m_{π}^{2}$ and $\int_{0}^{1} x (1 - x) d x = 1/6$ , $\int_{0}^{1} x^{2} (1 - x)^{2} d x = 1/30$ ,

L = 1 + \frac{t}{12} - \frac{t ^{2}}{120} + O (t^{3}) .

At $q = 1$ MeV, $t^{2} /120 = 2.3 \times 1 0^{- 11}$ , and the tape matches the two terms to $1 0^{- 13}$ .

A third representation is the dispersion integral: $L$ has an imaginary part above the two-pion threshold $μ = 2 m_{π}$ , and

L (q) = 1 + q^{2} \int_{2 m_{π}}^{\infty} \frac{μ ^{2} - 4 m _{π}^{2}}{μ ^{2} ( μ ^{2} + q ^{2} )} d μ .

The program maps $[2 m_{π}, \infty)$ to $[0, 1)$ with $μ = 2 m_{π} / (1 - u)$ and integrates by Simpson's rule on 2000 intervals. The integrand has a square-root zero at threshold, where Simpson's rule loses its order, and the result is off by $4.4 \times 1 0^{- 7}$ at $q = 200$ MeV.

[3]lines 209–299

Q_REF = [1.0, 50.0, 100.0, 200.0, 300.0, 500.0, 700.0, 1000.0]
L_REF = [1.0000043733362086695, 1.0107925820719110433, 1.0415967652983725786, 1.1464228856495893426, 1.2792753525586836379, 1.5487579788258103122, 1.7848080102098173246, 2.0735239825621998704]
L_tape = zeros(len(Q_REF))
for j to len(Q_REF) {
    L_tape[j] = eval_at(L_root, Q_REF[j])
}
show L_tape

#| The dispersion integrand after the map mu = 2 m_pi / (1 - u).
def disp_integrand(u, qq) {
    thr = 2.0 * MPI
    mu = thr / (1.0 - u)
    dmu = thr / ((1.0 - u) * (1.0 - u))
    return sqrt(mu * mu - thr * thr) / (mu * mu * (mu * mu + qq * qq)) * dmu
}

# Simpson intervals (even) for the dispersion integral
N_SIMPSON = 2000

#| L(q) from the once-subtracted dispersion integral, by Simpson's rule on N_SIMPSON intervals.
def L_disp(qq) {
    b = 1.0 - 1.0e-9
    h = b / N_SIMPSON
    s = disp_integrand(0.0, qq) + disp_integrand(b, qq)
    for i from 1 to N_SIMPSON {
        wt = 2.0
        if i % 2 == 1 {
            wt = 4.0
        }
        s = s + wt * disp_integrand(i * h, qq)
    }
    return 1.0 + qq * qq * s * h / 3.0
}
L_d = L_disp(200.0)
show L_d
disp_err = abs(L_d - L_REF[3])
show disp_err
L1_minus_one = L_tape[0] - 1.0
show L1_minus_one
t_1mev = 1.0 / (MPI * MPI)
t_sq_term = t_1mev * t_1mev / 120.0
show t_sq_term

# the two moments of the Feynman parameter in the small-q expansion, by
# Simpson's rule: exact for x(1 - x), and off by h^4 f''''/180 = 1.3e-13 for x^2 (1 - x)^2
N_MOMENT = 1000
h_m = 1.0 / N_MOMENT
moment_1 = 0.0
moment_2 = 0.0
for i to N_MOMENT + 1 {
    xm = i * h_m
    wm = 2.0
    if i % 2 == 1 {
        wm = 4.0
    }
    if i == 0 or i == N_MOMENT {
        wm = 1.0
    }
    moment_1 = moment_1 + wm * xm * (1.0 - xm)
    moment_2 = moment_2 + wm * xm * xm * (1.0 - xm) * (1.0 - xm)
}
moment_1 = moment_1 * h_m / 3.0
moment_2 = moment_2 * h_m / 3.0
show moment_1
show moment_2

test "the tape matches the Feynman-parameter integral at all eight momenta to 1e-14" {
    for j to len(Q_REF) {
        assert_near(L_tape[j], L_REF[j], 1.0e-14)
    }
}
test "at q = 1 MeV it is 1 + t/12 - t^2/120, t = q^2/m_pi^2, to 1e-13" {
    t1 = 1.0 / (MPI * MPI)
    assert_near(L_tape[0], 1.0 + t1 / 12.0 - t1 * t1 / 120.0, 1.0e-13)
}
test "the dispersion integral by Simpson's rule gives L(200) to 1e-6" {
    assert_near(L_d, 1.1464228856495893426, 1.0e-6)
}
test "the dispersion integral is off by 4.4e-7 at q = 200 MeV" {
    assert_near(disp_err, 4.42e-7, 0.01e-7)
}
test "at q = 1 MeV, L - 1 is 4.4e-6 and t^2/120 is 2.3e-11" {
    # to the tape's agreement with the integral, 1e-14
    assert_near(L1_minus_one, 1.0000043733362086695 - 1.0, 1.0e-14)
    # Python: (1 / 138.039**2)**2 / 120
    assert_near(t_sq_term, 2.2951524411056135e-11, 1.0e-24)
}
test "the integrals of x(1 - x) and x^2 (1 - x)^2 over [0, 1] are 1/6 and 1/30" {
    assert_near(moment_1, 1.0 / 6.0, 1.0e-14)
    assert_near(moment_2, 1.0 / 30.0, 2.0e-13)
}

Q_REF = [1.0, 50.0, 100.0, 200.0, 300.0, 500.0, 700.0, 1000.0]
L_REF = [1.0000043733362086695, 1.0107925820719110433, 1.0415967652983725786, 1.1464228856495893426, 1.2792753525586836379, 1.5487579788258103122, 1.7848080102098173246, 2.0735239825621998704]
L_tape = zeros(len(Q_REF))
for j to len(Q_REF) {
    L_tape[j] = eval_at(L_root, Q_REF[j])
}
show L_tape

#| The dispersion integrand after the map mu = 2 m_pi / (1 - u).
def disp_integrand(u, qq) {
    thr = 2.0 * MPI
    mu = thr / (1.0 - u)
    dmu = thr / ((1.0 - u) * (1.0 - u))
    return sqrt(mu * mu - thr * thr) / (mu * mu * (mu * mu + qq * qq)) * dmu
}

# Simpson intervals (even) for the dispersion integral
N_SIMPSON = 2000

#| L(q) from the once-subtracted dispersion integral, by Simpson's rule on N_SIMPSON intervals.
def L_disp(qq) {
    b = 1.0 - 1.0e-9
    h = b / N_SIMPSON
    s = disp_integrand(0.0, qq) + disp_integrand(b, qq)
    for i from 1 to N_SIMPSON {
        wt = 2.0
        if i % 2 == 1 {
            wt = 4.0
        }
        s = s + wt * disp_integrand(i * h, qq)
    }
    return 1.0 + qq * qq * s * h / 3.0
}
L_d = L_disp(200.0)
show L_d
disp_err = abs(L_d - L_REF[3])
show disp_err
L1_minus_one = L_tape[0] - 1.0
show L1_minus_one
t_1mev = 1.0 / (MPI * MPI)
t_sq_term = t_1mev * t_1mev / 120.0
show t_sq_term

# the two moments of the Feynman parameter in the small-q expansion, by
# Simpson's rule: exact for x(1 - x), and off by h^4 f''''/180 = 1.3e-13 for x^2 (1 - x)^2
N_MOMENT = 1000
h_m = 1.0 / N_MOMENT
moment_1 = 0.0
moment_2 = 0.0
for i to N_MOMENT + 1 {
    xm = i * h_m
    wm = 2.0
    if i % 2 == 1 {
        wm = 4.0
    }
    if i == 0 or i == N_MOMENT {
        wm = 1.0
    }
    moment_1 = moment_1 + wm * xm * (1.0 - xm)
    moment_2 = moment_2 + wm * xm * xm * (1.0 - xm) * (1.0 - xm)
}
moment_1 = moment_1 * h_m / 3.0
moment_2 = moment_2 * h_m / 3.0
show moment_1
show moment_2

test "the tape matches the Feynman-parameter integral at all eight momenta to 1e-14" {
    for j to len(Q_REF) {
        assert_near(L_tape[j], L_REF[j], 1.0e-14)
    }
}
test "at q = 1 MeV it is 1 + t/12 - t^2/120, t = q^2/m_pi^2, to 1e-13" {
    t1 = 1.0 / (MPI * MPI)
    assert_near(L_tape[0], 1.0 + t1 / 12.0 - t1 * t1 / 120.0, 1.0e-13)
}
test "the dispersion integral by Simpson's rule gives L(200) to 1e-6" {
    assert_near(L_d, 1.1464228856495893426, 1.0e-6)
}
test "the dispersion integral is off by 4.4e-7 at q = 200 MeV" {
    assert_near(disp_err, 4.42e-7, 0.01e-7)
}
test "at q = 1 MeV, L - 1 is 4.4e-6 and t^2/120 is 2.3e-11" {
    # to the tape's agreement with the integral, 1e-14
    assert_near(L1_minus_one, 1.0000043733362086695 - 1.0, 1.0e-14)
    # Python: (1 / 138.039**2)**2 / 120
    assert_near(t_sq_term, 2.2951524411056135e-11, 1.0e-24)
}
test "the integrals of x(1 - x) and x^2 (1 - x)^2 over [0, 1] are 1/6 and 1/30" {
    assert_near(moment_1, 1.0 / 6.0, 1.0e-14)
    assert_near(moment_2, 1.0 / 30.0, 2.0e-13)
}

The derivative

The derivative tape, evaluated at three momenta, matches mpmath's numerical derivative of the Feynman-parameter integral to $2 \times 1 0^{- 17}$ per MeV, rounding in a number of about 0.001. A central difference on the tape with a step of 1 MeV is off by up to $5.3 \times 1 0^{- 9}$ , which is the $h^{2} L^{'''} /6$ error of that step; the symbolic derivative has no step to choose.

[4]lines 311–368

Q_D = [50.0, 200.0, 500.0]
DL_REF = [0.00042617365415873921187, 0.0012407510640255769467, 0.0012762845936279545397]
dL_tape = zeros(len(Q_D))
dL_fd = zeros(len(Q_D))
H_FD = 1.0
for j to len(Q_D) {
    dL_tape[j] = eval_at(dL_root, Q_D[j])
    dL_fd[j] = (eval_at(L_root, Q_D[j] + H_FD) - eval_at(L_root, Q_D[j] - H_FD)) / (2.0 * H_FD)
}
show dL_tape
show dL_fd
dL_tape_err = 0.0
dL_fd_err = 0.0
for j to len(Q_D) {
    dL_tape_err = max(dL_tape_err, abs(dL_tape[j] - DL_REF[j]))
    dL_fd_err = max(dL_fd_err, abs(dL_fd[j] - DL_REF[j]))
}
show dL_tape_err
show dL_fd_err

test "the symbolic derivative matches mpmath's derivative of the integral to 1e-15 per MeV" {
    for j to len(Q_D) {
        assert_near(dL_tape[j], DL_REF[j], 1.0e-15)
    }
}
test "the derivative tape is within 2e-17 per MeV of mpmath at all three momenta" {
    assert_lt(dL_tape_err, 2.0e-17)
}
test "the central difference is off by up to 5.3e-9" {
    assert_near(dL_fd_err, 5.29e-9, 0.01e-9)
}
test "a central difference with a 1 MeV step is off by 2e-10 to 6e-9" {
    for j to len(Q_D) {
        assert_lt(abs(dL_fd[j] - DL_REF[j]), 6.0e-9)
        assert_lt(2.0e-10, abs(dL_fd[j] - DL_REF[j]))
    }
}

# L(q) and dL/dq from the tapes, on a grid of momenta
N_CURVE = 120
qs = linspace(1.0, 1000.0, N_CURVE)
Lc = zeros(N_CURVE)
dLc = zeros(N_CURVE)
for j to N_CURVE {
    Lc[j] = eval_at(L_root, qs[j])
    dLc[j] = 1000.0 * eval_at(dL_root, qs[j])
}
figure("The two-pion loop function $L(q)$, $m_\\pi = 138.039$ MeV", "momentum transfer q (MeV)", "L(q)")
plot(qs, Lc, "L(q), tape")
scatter(Q_REF, L_REF, "Feynman-parameter integral (mpmath)")
plot(qs, dLc, "dL/dq, derivative tape (per GeV)")
yright("dL/dq, derivative tape (per GeV)", "dL/dq (1/GeV)")
style("dL/dq, derivative tape (per GeV)", "dash")
scatter(Q_D, DL_REF * 1000.0, "mpmath derivative (per GeV)")
yright("mpmath derivative (per GeV)", "dL/dq (1/GeV)")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

Q_D = [50.0, 200.0, 500.0]
DL_REF = [0.00042617365415873921187, 0.0012407510640255769467, 0.0012762845936279545397]
dL_tape = zeros(len(Q_D))
dL_fd = zeros(len(Q_D))
H_FD = 1.0
for j to len(Q_D) {
    dL_tape[j] = eval_at(dL_root, Q_D[j])
    dL_fd[j] = (eval_at(L_root, Q_D[j] + H_FD) - eval_at(L_root, Q_D[j] - H_FD)) / (2.0 * H_FD)
}
show dL_tape
show dL_fd
dL_tape_err = 0.0
dL_fd_err = 0.0
for j to len(Q_D) {
    dL_tape_err = max(dL_tape_err, abs(dL_tape[j] - DL_REF[j]))
    dL_fd_err = max(dL_fd_err, abs(dL_fd[j] - DL_REF[j]))
}
show dL_tape_err
show dL_fd_err

test "the symbolic derivative matches mpmath's derivative of the integral to 1e-15 per MeV" {
    for j to len(Q_D) {
        assert_near(dL_tape[j], DL_REF[j], 1.0e-15)
    }
}
test "the derivative tape is within 2e-17 per MeV of mpmath at all three momenta" {
    assert_lt(dL_tape_err, 2.0e-17)
}
test "the central difference is off by up to 5.3e-9" {
    assert_near(dL_fd_err, 5.29e-9, 0.01e-9)
}
test "a central difference with a 1 MeV step is off by 2e-10 to 6e-9" {
    for j to len(Q_D) {
        assert_lt(abs(dL_fd[j] - DL_REF[j]), 6.0e-9)
        assert_lt(2.0e-10, abs(dL_fd[j] - DL_REF[j]))
    }
}

# L(q) and dL/dq from the tapes, on a grid of momenta
N_CURVE = 120
qs = linspace(1.0, 1000.0, N_CURVE)
Lc = zeros(N_CURVE)
dLc = zeros(N_CURVE)
for j to N_CURVE {
    Lc[j] = eval_at(L_root, qs[j])
    dLc[j] = 1000.0 * eval_at(dL_root, qs[j])
}
figure("The two-pion loop function $L(q)$, $m_\\pi = 138.039$ MeV", "momentum transfer q (MeV)", "L(q)")
plot(qs, Lc, "L(q), tape")
scatter(Q_REF, L_REF, "Feynman-parameter integral (mpmath)")
plot(qs, dLc, "dL/dq, derivative tape (per GeV)")
yright("dL/dq, derivative tape (per GeV)", "dL/dq (1/GeV)")
style("dL/dq, derivative tape (per GeV)", "dash")
scatter(Q_D, DL_REF * 1000.0, "mpmath derivative (per GeV)")
yright("mpmath derivative (per GeV)", "dL/dq (1/GeV)")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

The tape of L(q)

Values against the Feynman-parameter integral

The derivative

The tape of $L (q)$