Symbolic Mathematics
Dear reader, this chapter is a love letter to Lisp’s roots in symbolic computation. My interest in symbolic math software started in the 1980s when I had the Reduce symbolic math package running on my Xerox 1108 Lisp Machine. We build a small but complete symbolic mathematics library in Hy covering polynomial data structures, symbolic differentiation, and symbolic integration. I ported this example from the Common Lisp version in my book Loving Common Lisp, The Savvy Programmer’s Secret Weapon.
The beauty of this example is that Hy’s Lisp syntax makes the mathematical transformations read almost like the textbook rules they implement. The power rule for differentiation, for instance, is a single if expression. The fundamental theorem of calculus is three lines. This is exactly the kind of problem where Lisp shines.
Project Overview
The library is split into three modules that build on each other:
1 symbolic-math/
2 ├── data.hy ← Core classes: variables, constants, terms, polynomials, int\
3 egrals
4 ├── differentiation.hy ← Symbolic differentiation (power rule, sum rule)
5 ├── integration.hy ← Symbolic integration (reverse power rule, FTC)
6 ├── main.hy ← Combined smoke test runner
7 └── pyproject.toml ← Dependencies (just hy)
To run the complete test suite:
1 uv sync
2 uv run hy main.hy
The Data Layer: data.hy
The foundation of any symbolic math system is its data representation. We need types for variables, constants, terms (monomials), polynomials, and integrals. In Common Lisp these were defstruct definitions; in Hy we use defclass.
Symbolic Variables
A variable has a name and a domain (real, complex, or integer):
1 (defclass SymVariable []
2 "A symbolic variable such as x, y, or t."
3 (defn __init__ [self name [domain "real"]]
4 (setv self.name name)
5 (setv self.domain domain))
6
7 (defn __eq__ [self other]
8 (and (isinstance other SymVariable)
9 (= self.name other.name)
10 (= self.domain other.domain))))
We also define a constructor function and a predicate following Common Lisp naming conventions:
1 (defn make-variable [name [domain "real"]]
2 "Create a symbolic variable with name and optional domain."
3 (assert (in domain ["real" "complex" "integer"]))
4 (SymVariable name domain))
5
6 (defn variable-p [obj]
7 "Return True if OBJ is a SymVariable."
8 (isinstance obj SymVariable))
The -p suffix on variable-p is a Lisp tradition for predicate functions (from “predicate”). And make-variable wraps the constructor with validation, just as the Common Lisp version used check-type assertions.
Symbolic Constants
Constants like π and e can be stored symbolically and resolved to numeric values when needed:
1 (defclass SymConstant []
2 "A named mathematical constant."
3 (defn __init__ [self name value]
4 (setv self.name name)
5 (setv self.value value))
6
7 (defn numeric-value [self]
8 (cond
9 (= self.value "pi") math.pi
10 (= self.value "e") (math.exp 1.0)
11 True self.value)))
12
13 (defn make-constant [name value]
14 (SymConstant name value))
This lets us write integral bounds like (make-constant "pi" "pi") and defer numeric evaluation until the moment we actually need a floating-point result.
Terms (Monomials)
A term represents a single monomial: coefficient × variable^exponent. For example, 3x² is the term with coefficient 3, variable x, and exponent 2:
1 (defclass SymTerm []
2 "A single monomial: COEFFICIENT * VARIABLE ^ EXPONENT."
3 (defn __init__ [self coefficient variable exponent]
4 (setv self.coefficient coefficient)
5 (setv self.variable variable)
6 (setv self.exponent exponent))
7
8 (defn negate [self]
9 (SymTerm (- self.coefficient) self.variable self.exponent))
10
11 (defn scale [self scalar]
12 (SymTerm (* scalar self.coefficient) self.variable self.exponent))
13
14 (defn to-string [self]
15 (cond
16 (= self.exponent 0) (str self.coefficient)
17 (= self.exponent 1) f"{self.coefficient}{self.variable.name}"
18 True f"{self.coefficient}{self.variable.name}^{self.exponent}")))
Note how the to-string method handles three cases: constant terms (exponent 0), linear terms (exponent 1, no caret), and general terms. This mirrors exactly the Common Lisp term->string function.
Polynomials
A polynomial is a list of terms in one variable, kept sorted in descending exponent order. The key design decision, following the Common Lisp original, is that the make-polynomial constructor automatically combines like terms and drops zero-coefficient terms:
1 (defn %sort-and-combine-terms [terms]
2 "Internal: combine like-exponent terms, drop zero coefficients, sort descending."
3 (setv table {})
4 (for [term terms]
5 (setv e term.exponent)
6 (setv (get table e) (+ (.get table e 0) term.coefficient)))
7 (setv result [])
8 (for [[exp coeff] (.items table)]
9 (when (not (math.isclose coeff 0.0 :abs-tol 1e-12))
10 (.append result #(exp coeff))))
11 (sorted result :key (fn [x] (get x 0)) :reverse True))
This function uses a dictionary as an accumulator (keyed by exponent) to merge coefficients. The % prefix is a Lisp convention for internal/private helper functions.
With the normalizer in place, polynomial arithmetic becomes straightforward:
1 (defn polynomial-add [p q]
2 "Return P + Q as a new polynomial."
3 (make-polynomial p.variable (+ p.terms q.terms) :domain p.domain))
4
5 (defn polynomial-negate [poly]
6 (make-polynomial poly.variable
7 (lfor term poly.terms (term-negate term))
8 :domain poly.domain))
9
10 (defn polynomial-subtract [p q]
11 (polynomial-add p (polynomial-negate q)))
12
13 (defn polynomial-scale [poly scalar]
14 (make-polynomial poly.variable
15 (lfor term poly.terms (term-scale term scalar))
16 :domain poly.domain))
Addition simply concatenates the term lists and lets the normalizer handle merging. Subtraction negates then adds. This immutable, value oriented approach means we never mutate a polynomial and every operation returns a fresh one.
Evaluation substitutes a numeric value for the variable:
1 (defn polynomial-evaluate [poly value]
2 "Evaluate POLY at value."
3 (sum (lfor term poly.terms
4 (* term.coefficient (** value term.exponent)))))
Integrals
The integral structure wraps a polynomial with optional bounds, enabling both indefinite and definite integral representations:
1 (defn make-integral [integrand variable [lower None] [upper None]]
2 "Create a symbolic integral."
3 (assert (or (and (is lower None) (is upper None))
4 (and (not (is lower None)) (not (is upper None))))
5 "Either both lower and upper must be provided, or neither.")
6 (SymIntegral integrand variable lower upper))
The string representation uses Unicode for a clean mathematical display:
1 (defn integral->string [integral]
2 (setv var (str (variable-name integral.variable)))
3 (setv body (polynomial->string integral.integrand))
4 (setv bounds
5 (if (integral-definite-p integral)
6 f"[{(%bound->string integral.lower)},{(%bound->string integral.upper)}]"
7 ""))
8 f"∫{bounds}({body}) d{var}")
Data Layer Smoke Test
Running uv run hy data.hy exercises the data layer:
1 === Symbolic Math Data Layer Smoke Test ===
2
3 p : 3x^2 + -1x + 5
4 q : 1x + 2
5 p + q : 3x^2 + 7
6 p - q : 3x^2 + -2x + 3
7 2 * p : 6x^2 + -2x + 10
8 degree(p) : 2
9 p(0) : 5
10 p(1) : 7
11 p(2) : 15
12 indef : ∫(3x^2 + -1x + 5) dx
13 def : ∫[0,1](3x^2 + -1x + 5) dx
14 pi-int : ∫[0,pi](3x^2 + -1x + 5) dx
15
16 ===========================================
Notice how p + q correctly combined the x terms: -1x from p and 1x from q cancel out, leaving 3x^2 + 7.
Symbolic Differentiation: differentiation.hy
Calculus students learn three rules that suffice for differentiating any polynomial:
- Power rule: d/dx(c · xn) = n · c · x(n-1)
- Sum rule: d/dx(p + q) = p’ + q’
- Constant rule: d/dx(c) = 0
In Hy, the power rule is a single function:
1 (defn differentiate-term [term]
2 "Apply the power rule to a single monomial TERM."
3 (setv c term.coefficient)
4 (setv v term.variable)
5 (setv n term.exponent)
6 (if (= n 0)
7 None
8 (make-term (* n c) v (- n 1))))
When the exponent is zero (a constant term), the derivative is zero—we return None and filter it out at the polynomial level. Otherwise, we multiply the coefficient by the exponent and reduce the exponent by one. That’s the entire power rule.
Polynomial differentiation maps this over the terms and lets the constructor normalize:
1 (defn differentiate [poly]
2 "Differentiate polynomial POLY with respect to its variable."
3 (setv diff-terms [])
4 (for [term poly.terms]
5 (setv dt (differentiate-term term))
6 (when (not (is dt None))
7 (.append diff-terms dt)))
8 (if (not diff-terms)
9 (zero-polynomial poly.variable)
10 (make-polynomial poly.variable diff-terms :domain poly.domain)))
Higher-order derivatives use natural recursion:
1 (defn differentiate-n [poly n]
2 "Return the N-th derivative of polynomial POLY."
3 (if (= n 0)
4 poly
5 (differentiate-n (differentiate poly) (- n 1))))
We also provide convenience functions for numerical evaluation:
1 (defn gradient-at [poly value]
2 "Return the numerical value of the derivative of POLY at VALUE."
3 (polynomial-evaluate (differentiate poly) value))
4
5 (defn critical-point-p [poly value [tolerance 1e-9]]
6 "Return True if VALUE is (numerically) a critical point of POLY."
7 (< (abs (gradient-at poly value)) tolerance))
Differentiation Smoke Test
Running uv run hy differentiation.hy:
1 === Differentiation Smoke Test ===
2
3 p : 3x^2 + -1x + 5
4 p' : 6x + -1
5 p'' : 6
6 p''' : 0
7
8 q : 1x^4 + -2x^3 + 1x
9 q' : 4x^3 + -6x^2 + 1
10 q'' (via n=2): 12x^2 + -12x
11
12 gradient-at(p, 0) : -1 (expected -1)
13 gradient-at(p, 1) : 5 (expected 5)
14 critical-point-p(p,1/6): True (expected True)
15 critical-point-p(p,0) : False (expected False)
16
17 ==================================
The critical point test confirms that p’(1/6) = 0, since p’(x) = 6x - 1 and 6(1/6) - 1 = 0.
Symbolic Integration: integration.hy
Integration is the reverse of differentiation. For polynomials, the key rule is the reverse power rule:
∫ c · xn dx = (c / (n+1)) · x(n+1) + C
Notice the division by (n+1). In Common Lisp, dividing two integers produces a rational number automatically. Python’s integers don’t do this—3 / 2 gives 1.5, not 3/2. We use Python’s fractions.Fraction class to preserve exact rational coefficients:
1 (import fractions [Fraction])
2
3 (defn integrate-term [term]
4 "Apply the reverse power rule to a single monomial TERM."
5 (setv c term.coefficient)
6 (setv v term.variable)
7 (setv n term.exponent)
8 (setv new-exp (+ n 1))
9 (setv new-coef (if (isinstance c #(int Fraction))
10 (Fraction c new-exp)
11 (/ c new-exp)))
12 (make-term new-coef v new-exp))
When the coefficient is an integer or already a Fraction, we produce a Fraction result—keeping the representation exact. This means integrating 3x² gives 1x³ (since Fraction(3, 3) simplifies to 1), and integrating -1x gives -1/2x² (since Fraction(-1, 2) stays as-is).
The polynomial antiderivative maps the term integrator:
1 (defn integrate [poly]
2 "Return the antiderivative of polynomial POLY (without +C)."
3 (setv new-terms (lfor term poly.terms (integrate-term term)))
4 (if (not new-terms)
5 (zero-polynomial poly.variable)
6 (make-polynomial poly.variable new-terms :domain poly.domain)))
The Fundamental Theorem of Calculus
For definite integrals, we apply the fundamental theorem: ∫[a,b] f(x) dx = F(b) - F(a), where F is the antiderivative.
1 (defn evaluate-definite [poly lower upper]
2 "Evaluate the definite integral ∫[lower,upper] POLY dx numerically."
3 (setv F (integrate poly))
4 (setv a (%resolve-bound lower))
5 (setv b (%resolve-bound upper))
6 (- (polynomial-evaluate F b)
7 (polynomial-evaluate F a)))
The %resolve-bound helper converts bounds—which may be plain numbers or symbolic constants like π—to floating-point values for evaluation.
Integration Smoke Test
Running uv run hy integration.hy:
1 === Integration Smoke Test ===
2
3 p : 3x^2 + -1x + 5
4 ∫p dx : 1x^3 + -1/2x^2 + 5x
5 ∫∫p dx dx : 1/4x^4 + -1/6x^3 + 5/2x^2
6
7 q : 6x + 2
8 ∫q dx : 3x^2 + 2x
9
10 ∫₀¹ p dx : 5.5 (expected 5.5)
11 ∫₀¹ q dx : 5.0 (expected 5.0)
12 ∫₀^π p dx : 41.7794377477041
13
14 indef shell : ∫(3x^2 + -1x + 5) dx
15 def [0,1] : ∫[0,1](3x^2 + -1x + 5) dx
16 def [0,pi] : ∫[0,pi](3x^2 + -1x + 5) dx
17
18 ==============================
Notice the exact rational coefficients: -1/2, 1/4, -1/6, 5/2. These come from Fraction arithmetic and match what Common Lisp produces natively. When we evaluate the definite integral ∫₀¹ p dx, the antiderivative F(x) = x³ - (1/2)x² + 5x is evaluated at x=1 and x=0, giving F(1) - F(0) = 1 - 0.5 + 5 - 0 = 5.5.
Design Notes: From Common Lisp to Hy
Porting this library revealed several interesting contrasts between Common Lisp and Hy/Python:
Structs vs. Classes. Common Lisp’s defstruct generates constructors, accessors, and predicates automatically. In Hy we define classes with explicit __init__ and __eq__ methods, plus standalone accessor functions like term-coefficient for API compatibility with the original.
Exact Rationals. Common Lisp keeps rationals exact by default—(/ 3 2) returns 3/2, not 1.5. Python’s integers don’t do this, so we import fractions.Fraction to preserve exact coefficients during integration. This is the one area where Common Lisp is genuinely more convenient for symbolic math.
Immutability by Convention. Both versions follow the same discipline: every arithmetic operation returns a new polynomial; nothing is ever mutated in place. In Common Lisp this is natural with defstruct; in Python/Hy we simply choose not to modify attributes after construction.
Naming Conventions. We kept the Lisp-style names—make-polynomial, variable-p, polynomial->string—which Hy auto-translates to Python’s make_polynomial, variable_p, polynomial_to_string. The -p suffix for predicates and -> for conversions are idiomatic Lisp that read naturally in Hy.
Running the Full Test Suite
The file main.hy ties everything together:
1 (import data)
2 (import differentiation)
3 (import integration)
4
5 (defn main []
6 (print "=========================================================")
7 (print "=== RUNNING ALL SYMBOLIC MATH SMOKE TESTS IN HY ===")
8 (print "=========================================================")
9 (data.run-smoke-test)
10 (differentiation.run-smoke-test)
11 (integration.run-smoke-test))
12
13 (when (= __name__ "__main__")
14 (main))
Run it with:
1 uv run hy main.hy
Summary
We built a symbolic mathematics library from scratch in Hy, porting the data structures and algorithms from a Common Lisp implementation. The library demonstrates:
- Lisp-style data modeling with classes and constructor/predicate functions
- Immutable arithmetic where every operation returns fresh structures
-
Exact rational arithmetic using Python’s
fractions.Fraction - Recursive algorithms for higher-order derivatives and iterated integrals
- The fundamental theorem of calculus implemented in three lines of code
This is the kind of problem that Lisp was born to solve. The symbolic transformations—differentiating a term, integrating a polynomial, evaluating a definite integral—read almost identically to the mathematical rules they implement. Hy gives us this expressiveness while staying inside Python’s ecosystem.
If you enjoy symbolic math software, you can try extending this example code:
-
Polynomial multiplication — Implement
polynomial-multiplythat distributes each term of one polynomial across the terms of another, then normalizes. -
Pretty printing — Modify
term->stringto suppress the coefficient when it is 1 or -1 (printingx^2instead of1x^2and-xinstead of-1x). -
Chain rule — Extend differentiation to handle compositions like
(3x + 1)^4by representing composite expressions as a new data type. - Numeric integration — Implement Simpson’s rule as an alternative to the exact antiderivative approach, and compare accuracy for different step sizes.