Notebooks · lib/cas/check/pion_loop

The pion-loop function's derivative by lib/cas

posts/2026-09-25-the-pion-loop-function-by-computer-algebra.viv built the two-pion loop function $L (q) = \frac{w}{q} ln \frac{w + q}{2 m _{π}}$ , $w = 4 m_{π}^{2} + q^{2}$ , as an expression tape of 11 entries and differentiated it into a tape of 36, because nothing on the tape is shared or simplified. Here the same function is read into lib/cas from its formula, cas_diff differentiates it, and both are held to that post's references: a 30-digit Feynman-parameter integral for $L$ and mpmath's derivative of it for $d L / d q$ (posts/refs/the-pion-loop-function-by-computer-algebra.py), the same numbers to the same tolerances.

Three more routes to the derivative at 200 MeV agree with the CAS's: the differentiate phrase (forward mode), the same phrase using the symbolic derivative, and grad (reverse mode, lang/viv/src/grad.rs) of a def that computes $L$ . grad.rs and cas_diff take the derivatives of sqrt and log from one table.

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 22–86

import cas.cas as cas_module
from cas.calculus import cas_dag_size, cas_diff, cas_subst
from cas.parse import cas
from cas.store import cas_sym
from cas.eval import cas_value_at
from cas.show import cas_tex, cas_text

MPI = 138.039

L_TEXT = "sqrt(4*m^2 + q^2)/q * log((sqrt(4*m^2 + q^2) + q)/(2*m))"
q = cas_sym("q")
L = cas_subst(cas(L_TEXT), cas_sym("m"), cas("138.039"))
dL = cas_diff(L, q)
n_L = cas_dag_size(L)
n_dL = cas_dag_size(dL)
show n_L, n_dL
print("dL/dq = " + cas_text(dL))
print("LaTeX: " + cas_tex(dL))

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]
Q_D = [50.0, 200.0, 500.0]
DL_REF = [0.00042617365415873921187, 0.0012407510640255769467, 0.0012762845936279545397]

L_cas = zeros(len(Q_REF))
for j to len(Q_REF) {
    L_cas[j] = cas_value_at(L, q, Q_REF[j])
}
dL_cas = zeros(len(Q_D))
dL_err = 0.0
for j to len(Q_D) {
    dL_cas[j] = cas_value_at(dL, q, Q_D[j])
    dL_err = max(dL_err, abs(dL_cas[j] - DL_REF[j]))
}
show L_cas, dL_cas, dL_err

#| L(q) as a program computes it, for grad.
def loop_l(qq) {
    w = sqrt(4.0 * MPI * MPI + qq * qq)
    return w / qq * log((w + qq) / (2.0 * MPI))
}
d_forward = differentiate sqrt(4*MPI^2 + x^2)/x * log((sqrt(4*MPI^2 + x^2) + x)/(2*MPI)) at 200
d_symbolic = differentiate sqrt(4*MPI^2 + x^2)/x * log((sqrt(4*MPI^2 + x^2) + x)/(2*MPI)) at 200 using the symbolic derivative
d_reverse = grad(loop_l, 200.0)
show d_forward, d_symbolic, d_reverse

test "L from the CAS matches the Feynman-parameter integral at all eight momenta to 1e-14" {
    for j to len(Q_REF) {
        assert_near(L_cas[j], L_REF[j], 1.0e-14)
    }
}
test "dL/dq from cas_diff matches mpmath's derivative to 2e-17 per MeV, as the post's tape does" {
    assert_lt(dL_err, 2.0e-17)
}
test "the derivative's DAG is smaller than the post's 36-entry tape" {
    assert_lt(n_dL, 36)
}
test "forward mode, the symbolic phrase and reverse mode agree with cas_diff at 200 MeV" {
    assert_near(d_forward, dL_cas[1], 1.0e-17)
    assert_near(d_symbolic, dL_cas[1], 1.0e-17)
    assert_near(d_reverse, dL_cas[1], 1.0e-17)
}

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

import cas.cas as cas_module
from cas.calculus import cas_dag_size, cas_diff, cas_subst
from cas.parse import cas
from cas.store import cas_sym
from cas.eval import cas_value_at
from cas.show import cas_tex, cas_text

MPI = 138.039

L_TEXT = "sqrt(4*m^2 + q^2)/q * log((sqrt(4*m^2 + q^2) + q)/(2*m))"
q = cas_sym("q")
L = cas_subst(cas(L_TEXT), cas_sym("m"), cas("138.039"))
dL = cas_diff(L, q)
n_L = cas_dag_size(L)
n_dL = cas_dag_size(dL)
show n_L, n_dL
print("dL/dq = " + cas_text(dL))
print("LaTeX: " + cas_tex(dL))

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]
Q_D = [50.0, 200.0, 500.0]
DL_REF = [0.00042617365415873921187, 0.0012407510640255769467, 0.0012762845936279545397]

L_cas = zeros(len(Q_REF))
for j to len(Q_REF) {
    L_cas[j] = cas_value_at(L, q, Q_REF[j])
}
dL_cas = zeros(len(Q_D))
dL_err = 0.0
for j to len(Q_D) {
    dL_cas[j] = cas_value_at(dL, q, Q_D[j])
    dL_err = max(dL_err, abs(dL_cas[j] - DL_REF[j]))
}
show L_cas, dL_cas, dL_err

#| L(q) as a program computes it, for grad.
def loop_l(qq) {
    w = sqrt(4.0 * MPI * MPI + qq * qq)
    return w / qq * log((w + qq) / (2.0 * MPI))
}
d_forward = differentiate sqrt(4*MPI^2 + x^2)/x * log((sqrt(4*MPI^2 + x^2) + x)/(2*MPI)) at 200
d_symbolic = differentiate sqrt(4*MPI^2 + x^2)/x * log((sqrt(4*MPI^2 + x^2) + x)/(2*MPI)) at 200 using the symbolic derivative
d_reverse = grad(loop_l, 200.0)
show d_forward, d_symbolic, d_reverse

test "L from the CAS matches the Feynman-parameter integral at all eight momenta to 1e-14" {
    for j to len(Q_REF) {
        assert_near(L_cas[j], L_REF[j], 1.0e-14)
    }
}
test "dL/dq from cas_diff matches mpmath's derivative to 2e-17 per MeV, as the post's tape does" {
    assert_lt(dL_err, 2.0e-17)
}
test "the derivative's DAG is smaller than the post's 36-entry tape" {
    assert_lt(n_dL, 36)
}
test "forward mode, the symbolic phrase and reverse mode agree with cas_diff at 200 MeV" {
    assert_near(d_forward, dL_cas[1], 1.0e-17)
    assert_near(d_symbolic, dL_cas[1], 1.0e-17)
    assert_near(d_reverse, dL_cas[1], 1.0e-17)
}

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