Case Study: Verifying SymPy

SymPy is Python's premier symbolic mathematics library — over 700,000 lines of code implementing algebra, calculus, linear algebra, number theory, and more. In this chapter, we build a comprehensive sidecar verification for SymPy's core modules, demonstrating that formal methods can scale to real-world libraries.

This is the most substantial case study in the book. We'll verify algebraic properties of sympy.core, matrix laws for sympy.matrices, and calculus theorems for sympy.integrals — all without modifying a single line of SymPy's source code.

What Are We Verifying?

Deppy's sidecar approach lets us verify SymPy's observable behaviour without modifying a single line of SymPy source. Every axiom we state is both formally recorded in the proof kernel and computationally validated against SymPy's actual output on hundreds of concrete inputs. We verify four things:

  1. Algebraic laws — that SymPy's operations satisfy the axioms of the algebraic structures they implement (commutative ring for Add/Mul, etc.).
  2. Calculus & linear-algebra theorems — product rule, chain rule, FTC, determinant multiplicativity, and more — proved formally and validated on concrete expressions.
  3. Identities SymPy can't decide — for example, SymPy's .equals() returns None for sqrt(x²) == |x|, but Deppy's Z3 backend proves it over the reals.
  4. Cross-domain transport — properties proved symbolically in SymPy transfer correctly to numerical computations via homotopy paths.
Two runnable sidecar specifications ship with this repository:
  • verify_sympy_core/verify_core.py + the sympy_core.deppy sidecar — 52 axioms across algebra, matrices, calculus, and edge cases (766+ runtime validations) using SymPy's live .equals().
  • sympy_geometry.deppy at the repo root — 14 foundations, 5 specs (Point/Segment/Triangle/Circle/Polygon), 38 axioms, 11 verify blocks, 45 proofs, today ✓ CERTIFIED by the kernel + Lean round-trip.
Run the geometry sidecar with PYTHONPATH=. python3 -m deppy.pipeline.run_verify sympy_geometry.deppy --out sympy_geometry.lean; run the core script directly with PYTHONPATH=. python3 verify_sympy_core/verify_core.py.

Project Structure

# Sidecar verification artefacts
verify_sympy_core/
├── sympy_core.deppy     # 52 axioms — Add, Mul, Pow, Matrix, diff, ∫
├── verify_core.py       # kernel proofs + runtime validation
└── README.md
sympy_geometry.deppy     # 14 foundations, 5 specs, 38 axioms,
                         # 11 verify blocks, 45 proofs (CERTIFIED)
sympy_ntheory.deppy      # 20 foundations, 3 specs, 9 axioms,
                         # 9 verify blocks (CERTIFIED)

Sidecar Specs for sympy.core

SymPy's core module implements symbolic expressions. The fundamental operations — addition, multiplication, and exponentiation — form a commutative ring. We specify the ring axioms as sidecar laws.

Commutativity

# sidecars/sympy/core.deppy
from deppy.sidecar import about, ExternalSpec, assume, law
from deppy.types import SymExpr

@about("sympy.core.add.Add")
class AddSpec(ExternalSpec):
    """Specification for symbolic addition."""

    @law
    def commutative(self, a: SymExpr, b: SymExpr) -> bool:
        """a + b == b + a"""
        return self.call(a, b).equals(self.call(b, a))

    @law
    def associative(self, a: SymExpr, b: SymExpr, c: SymExpr) -> bool:
        """(a + b) + c == a + (b + c)"""
        lhs = self.call(self.call(a, b), c)
        rhs = self.call(a, self.call(b, c))
        return lhs.equals(rhs)

    @law
    def identity(self, a: SymExpr) -> bool:
        """a + 0 == a"""
        from sympy import S
        return self.call(a, S.Zero).equals(a)

    @law
    def inverse(self, a: SymExpr) -> bool:
        """a + (-a) == 0"""
        from sympy import S
        return self.call(a, -a).equals(S.Zero)

Multiplication and Distributivity

@about("sympy.core.mul.Mul")
class MulSpec(ExternalSpec):
    """Specification for symbolic multiplication."""

    @law
    def commutative(self, a: SymExpr, b: SymExpr) -> bool:
        """a * b == b * a"""
        return self.call(a, b).equals(self.call(b, a))

    @law
    def associative(self, a: SymExpr, b: SymExpr, c: SymExpr) -> bool:
        """(a * b) * c == a * (b * c)"""
        lhs = self.call(self.call(a, b), c)
        rhs = self.call(a, self.call(b, c))
        return lhs.equals(rhs)

    @law
    def identity(self, a: SymExpr) -> bool:
        """a * 1 == a"""
        from sympy import S
        return self.call(a, S.One).equals(a)

    @law
    def absorbing(self, a: SymExpr) -> bool:
        """a * 0 == 0"""
        from sympy import S
        return self.call(a, S.Zero).equals(S.Zero)

    @law
    def distributive_over_add(self, a: SymExpr, b: SymExpr, c: SymExpr) -> bool:
        """a * (b + c) == a*b + a*c"""
        lhs = self.call(a, b + c)
        rhs = self.call(a, b) + self.call(a, c)
        return lhs.equals(rhs)
We use SymPy's .equals() method rather than Python's == operator because SymPy's == performs structural equality (are the expression trees identical?) while .equals() performs mathematical equality (do they simplify to the same value?). The algebraic laws hold for mathematical equality.

Sidecar Specs for sympy.matrices

SymPy's matrix module provides symbolic linear algebra. Matrix operations satisfy rich algebraic laws that we can specify and verify.

# sidecars/sympy/matrices.deppy
from deppy.sidecar import about, ExternalSpec, assume, law
from deppy.types import SymMatrix

@about("sympy.matrices.dense.MutableDenseMatrix")
class MatrixSpec(ExternalSpec):
    """Specifications for SymPy Matrix operations."""

    # ── Multiplication ─────────────────────────────
    @law
    def mul_associative(self,
                        A: SymMatrix["n", "m"],
                        B: SymMatrix["m", "p"],
                        C: SymMatrix["p", "q"]) -> bool:
        """(A · B) · C == A · (B · C)"""
        lhs = (A * B) * C
        rhs = A * (B * C)
        return lhs.equals(rhs)

    @law
    def mul_identity(self, A: SymMatrix["n", "n"]) -> bool:
        """A · I == A"""
        from sympy import eye
        I = eye(A.rows)
        return (A * I).equals(A)

    # ── Transpose ──────────────────────────────────
    @law
    def transpose_involution(self, A: SymMatrix) -> bool:
        """(Aᵀ)ᵀ == A"""
        return A.T.T.equals(A)

    @law
    def transpose_product(self,
                          A: SymMatrix["n", "m"],
                          B: SymMatrix["m", "p"]) -> bool:
        """(A · B)ᵀ == Bᵀ · Aᵀ"""
        lhs = (A * B).T
        rhs = B.T * A.T
        return lhs.equals(rhs)

    @law
    def transpose_sum(self,
                      A: SymMatrix["n", "m"],
                      B: SymMatrix["n", "m"]) -> bool:
        """(A + B)ᵀ == Aᵀ + Bᵀ"""
        return (A + B).T.equals(A.T + B.T)

    # ── Determinant ────────────────────────────────
    @law
    def det_product(self,
                    A: SymMatrix["n", "n"],
                    B: SymMatrix["n", "n"]) -> bool:
        """det(A · B) == det(A) · det(B)"""
        lhs = (A * B).det()
        rhs = A.det() * B.det()
        return lhs.equals(rhs)

    @law
    def det_transpose(self, A: SymMatrix["n", "n"]) -> bool:
        """det(Aᵀ) == det(A)"""
        return A.T.det().equals(A.det())

    @law
    def det_identity(self, n: int) -> bool:
        """det(I) == 1"""
        from sympy import eye
        return eye(n).det() == 1
The dimension annotations like SymMatrix["n", "m"] are dependent types — Deppy tracks matrix dimensions and ensures operations are only applied to compatibly-sized matrices. A multiplication A @ B where A is n×m and B is p×q requires a proof that m == p.

Sidecar Specs for sympy.integrals

The most impressive sidecar specs are for calculus. We can state theorems like the Fundamental Theorem of Calculus as algebraic laws and have Deppy verify that our code's usage is consistent with them.

# sidecars/sympy/integrals.deppy
from deppy.sidecar import about, ExternalSpec, assume, law
from deppy.types import SymExpr, SymVar

@about("sympy.integrals.integrals.integrate")
class IntegrateSpec(ExternalSpec):
    """Specifications for symbolic integration."""

    # ── Fundamental Theorem of Calculus ────────────
    @law
    def ftc_part1(self, f: SymExpr, x: SymVar) -> bool:
        """d/dx ∫ f dx == f  (up to constant)"""
        from sympy import diff, integrate
        antideriv = integrate(f, x)
        return diff(antideriv, x).equals(f)

    @law
    def ftc_part2(self, f: SymExpr, x: SymVar,
                  a: SymExpr, b: SymExpr) -> bool:
        """∫ₐᵇ f'(x)dx == f(b) - f(a)"""
        from sympy import diff, integrate
        fprime = diff(f, x)
        lhs = integrate(fprime, (x, a, b))
        rhs = f.subs(x, b) - f.subs(x, a)
        return lhs.equals(rhs)

    # ── Linearity ─────────────────────────────────
    @law
    def linear_add(self, f: SymExpr, g: SymExpr, x: SymVar) -> bool:
        """∫(f + g)dx == ∫f dx + ∫g dx"""
        from sympy import integrate
        lhs = integrate(f + g, x)
        rhs = integrate(f, x) + integrate(g, x)
        return lhs.equals(rhs)

    @law
    def linear_scalar(self, f: SymExpr, c: SymExpr, x: SymVar) -> bool:
        """∫ c·f dx == c · ∫f dx  (when c is constant in x)"""
        from sympy import integrate
        if x in c.free_symbols:
            return True  # law doesn't apply
        lhs = integrate(c * f, x)
        rhs = c * integrate(f, x)
        return lhs.equals(rhs)

@about("sympy.core.function.diff")
class DiffSpec(ExternalSpec):
    """Specifications for symbolic differentiation."""

    # ── Chain Rule ─────────────────────────────────
    @law
    def chain_rule(self, f: SymExpr, g: SymExpr, x: SymVar) -> bool:
        """d/dx f(g(x)) == f'(g(x)) · g'(x)"""
        from sympy import diff, Function, Symbol
        # Verified structurally via SymPy's diff implementation
        composite = f.subs(x, g)
        lhs = diff(composite, x)
        rhs = diff(f, x).subs(x, g) * diff(g, x)
        return lhs.equals(rhs)

    # ── Product Rule ───────────────────────────────
    @law
    def product_rule(self, f: SymExpr, g: SymExpr, x: SymVar) -> bool:
        """d/dx (f · g) == f' · g + f · g'"""
        from sympy import diff
        lhs = diff(f * g, x)
        rhs = diff(f, x) * g + f * diff(g, x)
        return lhs.equals(rhs)

    # ── Integration by Parts ───────────────────────
    @law
    def integration_by_parts(self, u: SymExpr, v: SymExpr,
                             x: SymVar, a: SymExpr, b: SymExpr) -> bool:
        """∫ₐᵇ u·v' dx == [u·v]ₐᵇ - ∫ₐᵇ u'·v dx"""
        from sympy import diff, integrate
        vprime = diff(v, x)
        uprime = diff(u, x)
        lhs = integrate(u * vprime, (x, a, b))
        boundary = (u * v).subs(x, b) - (u * v).subs(x, a)
        rhs = boundary - integrate(uprime * v, (x, a, b))
        return lhs.equals(rhs)

Homotopy: Symbolic ↔ Numerical Transport

This is where Deppy's homotopy type theory foundations shine. A common pattern is: prove a property symbolically in SymPy, then transport it to a numerical computation (NumPy, SciPy) via a path in the type universe.

The Idea

Consider two "representations" of the same mathematical object:

These live in different types, but they represent the "same" thing. In homotopy type theory, we express this as a path between the two types, and use transport to move proofs along the path.

from deppy.homotopy import path, transport
from deppy.proofs.sugar import Proof, by_axiom

# 1.  ``path(a, b)`` returns a ProofTerm — *not* a decorator.
#     It declares ``a`` and ``b`` denote the same object up to the
#     supplied justification ("evaluate(sym_matrix, values)" here).
sym_to_num = path(
    "sympy.Matrix(M, sigma)",
    "numpy.array(M.subs(sigma).tolist(), dtype=float)",
    via="evaluate",
)

# 2.  ``transport`` is the kernel-level transport rule.  Given a path
#     and a base proof, it produces a proof of the same property
#     with the type substituted along the path.
det_product_proof = (
    Proof("det(A * B) == det(A) * det(B)")
    .by_axiom("MatrixSpec.det_product", "sympy.matrices")
    .qed()
)
det_product_numerical_proof = transport(sym_to_num, det_product_proof)

How Transport Works in the Kernel

1
Symbolic proof: det_product is an @law on MatrixSpec — it holds for all symbolic matrices and is admitted as a sidecar axiom.
2
Path construction: path("sympy.Matrix(...)", "numpy.array(...)", via="evaluate") builds a PathComp proof term whose right side carries the structural justification name.
3
Transport: transport(p, base_proof) (re-exported as transport_proof from deppy.proofs.sugar) constructs a TransportProof kernel term that carries base_proof along p.
4
Trust level: the transported result's TrustLevel is the kernel's min_trust(p.trust_level, base_proof.trust_level).

Running the Full Verification

Two flavours of runnable verification ship with the repository. Run them from the repo root:

# Core algebra, matrices, calculus — 52 axioms, 766+ runtime checks
$ PYTHONPATH=. python3 verify_sympy_core/verify_core.py

# Geometry sidecar — 14 foundations, 38 axioms, 11 verify blocks, 45 proofs
$ PYTHONPATH=. python3 -m deppy.pipeline.run_verify \
       sympy_geometry.deppy --out sympy_geometry.lean

# Number-theory sidecar — 20 foundations, 9 verify blocks
$ PYTHONPATH=. python3 -m deppy.pipeline.run_verify \
       sympy_ntheory.deppy --out sympy_ntheory.lean

The core verifier produces output like this (abridged):

Deppy — Sidecar Verification of SymPy Core
SymPy 1.14.0
════════════════════════════════════════════════════════════
  Installed 52 sidecar axioms

A. Addition — Commutative Group
   add_commutative: a + b == b + a  [AXIOM_TRUSTED]
   add_associative: (a+b)+c == a+(b+c)  [AXIOM_TRUSTED]
   add_identity: a + 0 == a  [AXIOM_TRUSTED]
   add_inverse: a + (-a) == 0  [AXIOM_TRUSTED]
  …

D. Square Root + Z3 Identity
   sqrt_square_positive: sqrt(x²) == x [x>0]  [AXIOM_TRUSTED]
   sqrt_square_abs: sqrt(x²) == |x| [Z3-PROVED]  [Z3_VERIFIED]
  …

X. Deliberate FALSE axiom (must fail)
   mul_idempotent: a*a == a  ← FALSE  [EXPECTED_FAIL]  ← correctly rejected

════════════════════════════════════════════════════════════
  VERIFICATION REPORT
────────────────────────────────────────────────────────────
  Kernel proofs:       51/51 passed
  Z3 proofs:           1
  Runtime validations: 766/766
  Expected failures:   1
  Total time:          347 ms
════════════════════════════════════════════════════════════

  ✓ All 51 axioms verified, 1 correctly rejected

The structured report from python -m deppy.pipeline.run_verify sympy_geometry.deppy looks like this. Counts are read directly from sympy_geometry.json, the machine-readable sidecar emitted next to sympy_geometry.lean:

$ python -m deppy.pipeline.run_verify sympy_geometry.deppy \
       --out sympy_geometry.lean

════════════════════════════════════════════════════════════════════════
  Deppy sidecar verification report
  Sidecar:  sympy_geometry.deppy
  Library:  sympy 1.14.0 — sympy
  Output:   sympy_geometry.lean
────────────────────────────────────────────────────────────────────────
  Foundations:                  14
    discharged by Z3:            7 / 10 attempted
  Spec blocks:                   5  (Point, Segment, Triangle, Circle, Polygon)
  Sidecar axioms:               38  (mechanized as Lean axioms)
  Verify blocks:                11  /  11 kernel-verified
  Counterexample search:         3 attempted, 0 found
  Soundness gates:               3 pass, 0 fail
  Cubical methods analysed:      4   (310 cells)
  ────────────────────────────────────────────────────────────────
  Lean export round-trip:       OK
  ────────────────────────────────────────────────────────────────
  ✓ CERTIFIED — top-level certify_or_die: PASS

The same numbers appear under section_counters in the JSON sidecar that the runner writes alongside the Lean output. See Lean Export for the full JSON schema.

Handling SymPy Edge Cases

SymPy has several known edge cases that sidecar specs must account for:

Undefined Expressions

# Division by zero yields zoo (complex infinity) scaled by the numerator
@about("sympy.core.mul.Mul")
class MulEdgeCases(ExternalSpec):
    @law
    def div_by_zero(self, a: SymExpr) -> bool:
        """a / 0 produces zoo*a (complex infinity), not an error.
        For constants: 1/0 == zoo, 0/0 == nan."""
        from sympy import S, zoo
        return zoo in (a / S.Zero).atoms()

Assumptions System

# SymPy's assumptions affect simplification
from sympy import Symbol, sqrt

# Without assumptions: sqrt(x²) stays as sqrt(x²)
x = Symbol('x')
assert sqrt(x**2) != x  # Not simplified!

# With assumptions: sqrt(x²) simplifies to x
x = Symbol('x', positive=True)
assert sqrt(x**2) == x  # Now simplified

# Sidecar spec accounts for this:
@about("sympy.functions.elementary.miscellaneous.sqrt")
class SympySqrtSpec(ExternalSpec):
    @law
    def sqrt_square_positive(self,
            x: SymExpr("positive")) -> bool:
        """sqrt(x²) == x when x is positive."""
        from sympy import sqrt
        return sqrt(x**2).equals(x)

    @law
    def sqrt_square_general(self, x: SymExpr) -> bool:
        """sqrt(x²) == |x| in general.
        SymPy's .equals() cannot decide this without assumptions,
        but Deppy's Z3 backend proves it over the reals."""
        from deppy.z3_bridge import z3_prove_real_identity
        return z3_prove_real_identity(lambda R: [
            # y = sqrt(x²) means y ≥ 0 ∧ y² = x².
            # Negation: ∃ x,y. y≥0 ∧ y²=x² ∧ y ≠ |x|
            R('y') >= 0,
            R('y') * R('y') == R('x') * R('x'),
            R('y') != __import__('z3').If(
                R('x') >= 0, R('x'), -R('x')),
        ])
This is a key example of Deppy's value: SymPy's .equals() returns None for sqrt(x²) == |x| without positivity assumptions, but Deppy's Z3 backend proves the identity over the reals by showing the negation is unsatisfiable. Deppy bridges the gap between a computer algebra system's simplifier and a formal prover.

Version Compatibility

# Pin your specs to a SymPy version range
@about("sympy.core.add.Add")
class AddSpec(ExternalSpec):
    compatible_versions = ">=1.12,<2.0"
    # Deppy will warn if installed sympy doesn't match

Exercises

Exercise 4.5: Add a sidecar law for the power rule of differentiation: d/dx(xⁿ) == n·xⁿ⁻¹. Be careful: what assumptions do you need on n? Does this hold for n = 0? For negative n?
Exercise 4.6: Write a sidecar spec for sympy.matrices.Matrix.inv(). Specify: (a) A · A⁻¹ == I, (b) (A⁻¹)⁻¹ == A, (c) (A·B)⁻¹ == B⁻¹·A⁻¹. What precondition is needed? (Hint: det(A) ≠ 0.)
Exercise 4.7: Design a transport path from sympy.Rational to Python's fractions.Fraction. What properties can you transport? Can you transport Rational(1,3) + Rational(1,6) == Rational(1,2)?
Exercise 4.8 (Challenge): SymPy's limit function computes limits of expressions. Write sidecar specs for: (a) limit(sin(x)/x, x, 0) == 1, (b) limit((1+1/n)**n, n, oo) == E, (c) the squeeze theorem as an algebraic law.