Symbolic Mathematics in Common Lisp
Dear reader, in the early 1970s I earned a Bachelor Of Science degree in Physics from UC Santa Barbara. Although little of my work in the last 50 years has involved either Physics or pure mathematics, I have a long term interest in symbolic math systems, starting when I installed the Reduce system on my Xerox Lisp Machine in 1982. Please note that the material here is for my own Lisp hacking enjoyment.
Common Lisp has long been a natural home for symbolic computation. Its homoiconic nature, rich macro system, and first-class support for rational arithmetic make it an ideal language for building systems that manipulate mathematical expressions as data rather than reducing them immediately to floating-point numbers. In this chapter we build a small but complete symbolic mathematics library from scratch. We define data structures for variables, constants, monomials, polynomials, and integrals; implement symbolic differentiation using the power rule; and implement symbolic integration using the reverse power rule together with the Fundamental Theorem of Calculus. The result is a readable, testable, purely functional library that illuminates both the mathematics and the Lisp idioms involved.
The Data Layer
Every symbolic computation system needs a representation for the objects it manipulates. Before we can differentiate or integrate anything we must decide how to store variables, constants, terms, polynomials, and integrals. The choices made here cascade through the rest of the system, so it is worth studying them carefully.
Design Philosophy
The library is deliberately CLOS-free. Every type is defined with defstruct, which gives us lightweight value structs with automatic constructors, slot accessors, and copy functions. All arithmetic operations return fresh structs; no existing object is ever mutated. This immutable, functional style keeps the code easy to reason about and makes it straightforward to use results as inputs to further computation without defensive copying.
The entire data layer lives in the SYMBOLIC-MATH package and is loaded from a single file, data.lisp. The package exports every public name: constructors, predicates, accessors, and helpers so that the differentiation and integration layers can import only what they need with (:use #:cl #:symbolic-math).
Variables and Constants
The two most primitive objects in the system are symbolic variables (unknowns like x or t) and named constants (fixed values like π or e). The following listing defines both structs together with their constructors, predicates, and a numeric-evaluation helper for constants.
1 (defstruct (sym-variable (:conc-name variable-))
2 "A symbolic variable such as x, y, or t.
3 DOMAIN may be :real (default), :complex, or :integer."
4 (name nil :type symbol)
5 (domain :real :type keyword))
6
7 (defun make-variable (name &key (domain :real))
8 "Create a symbolic variable with NAME (a symbol) and optional DOMAIN.
9 DOMAIN is one of :real (default), :complex, or :integer.
10
11 Example:
12 (make-variable 'x) ; => sym-variable x ∈ ℝ
13 (make-variable 't :domain :real)"
14
15 (check-type name symbol)
16 (assert (member domain '(:real :complex :integer))
17 (domain) "DOMAIN must be :real, :complex, or :integer, got ~s" domain)
18 (make-sym-variable :name name :domain domain))
19
20 (defun variable-p (obj)
21 "Return T if OBJ is a sym-variable."
22 (sym-variable-p obj))
23
24 (defun variable= (a b)
25 "Return T if variables A and B represent the same symbolic variable.
26 Two variables are equal when their names and domains match."
27 (and (sym-variable-p a)
28 (sym-variable-p b)
29 (eq (variable-name a) (variable-name b))
30 (eq (variable-domain a) (variable-domain b))))
31
32 (defstruct (sym-constant (:conc-name constant-))
33 "A named mathematical constant.
34 VALUE may be a Common Lisp number, :pi, or :e (Euler's number)."
35 (name nil :type symbol)
36 (value nil)) ; number | :pi | :e
37
38 (defun make-constant (name value)
39 "Create a named constant with NAME (symbol) and VALUE.
40 VALUE may be a real number, :pi, or :e."
41 (check-type name symbol)
42 (assert (or (numberp value) (member value '(:pi :e)))
43 (value) "VALUE must be a number, :pi, or :e, got ~s" value)
44 (make-sym-constant :name name :value value))
45
46 (defun constant-p (obj)
47 "Return T if OBJ is a sym-constant."
48 (sym-constant-p obj))
49
50 (defun constant-numeric-value (c)
51 "Return the numeric (floating-point) value of constant C.
52 :pi → pi, :e → e, numbers are returned as-is."
53 (check-type c sym-constant)
54 (let ((v (constant-value c)))
55 (cond
56 ((eq v :pi) pi)
57 ((eq v :e) (exp 1.0d0))
58 (t v))))
sym-variable stores a Lisp symbol as the variable name and a keyword denoting its domain. The (:conc-name variable-) option tells defstruct to prefix all generated accessor names with variable-, so the name slot is reached via variable-name rather than the default sym-variable-name. The make-variable wrapper validates its arguments with check-type and assert before delegating to the raw struct constructor, ensuring that only well-formed objects enter the system.
sym-constant follows the same pattern but its value slot is more flexible: it may hold a Common Lisp real number, the keyword :pi, or the keyword :e. The constant-numeric-value function resolves any of these to a double-float on demand using a short cond dispatch. Storing the symbolic name alongside the numeric value means that bounds such as π can be displayed in human-readable form (as shown later in the integral display code) while still being convertible to a floating-point number when a definite integral must be evaluated numerically.
Monomials (Terms)
A monomial is a single product of a coefficient, a variable, and a non-negative integer power is the atom from which polynomials are assembled. The sym-term struct and its helpers are shown below.
1 (defstruct (sym-term (:conc-name term-))
2 "A single monomial: COEFFICIENT * VARIABLE ^ EXPONENT.
3 EXPONENT is a non-negative integer (0 gives a constant term)."
4 (coefficient 0 :type real)
5 (variable nil) ; sym-variable
6 (exponent 0 :type (integer 0 *)))
7
8 (defun make-term (coefficient variable exponent)
9 "Create a monomial: COEFFICIENT * VARIABLE ^ EXPONENT.
10 EXPONENT must be a non-negative integer.
11
12 Examples:
13 (make-term 3 x 2) ; 3x²
14 (make-term -1 x 1) ; -x
15 (make-term 5 x 0) ; constant 5"
16
17 (check-type coefficient real)
18 (check-type variable sym-variable)
19 (check-type exponent (integer 0 *))
20 (make-sym-term :coefficient coefficient
21 :variable variable
22 :exponent exponent))
23
24 (defun term-p (obj)
25 "Return T if OBJ is a sym-term."
26 (sym-term-p obj))
27
28 (defun term= (a b)
29 "Return T if terms A and B are structurally equal
30 (same variable, exponent, and coefficient)."
31 (and (sym-term-p a) (sym-term-p b)
32 (= (term-coefficient a) (term-coefficient b))
33 (variable= (term-variable a) (term-variable b))
34 (= (term-exponent a) (term-exponent b))))
35
36 (defun term-negate (term)
37 "Return a new term equal to -TERM."
38 (make-sym-term :coefficient (- (term-coefficient term))
39 :variable (term-variable term)
40 :exponent (term-exponent term)))
41
42 (defun term-scale (term scalar)
43 "Return a new term equal to SCALAR * TERM."
44 (check-type scalar real)
45 (make-sym-term :coefficient (* scalar (term-coefficient term))
46 :variable (term-variable term)
47 :exponent (term-exponent term)))
48
49 (defun term->string (term)
50 "Return a human-readable string representation of TERM.
51 Example: 3x^2, -x^1, 5 (for exponent 0)."
52 (let ((c (term-coefficient term))
53 (v (symbol-name (variable-name (term-variable term))))
54 (e (term-exponent term)))
55 (cond
56 ((zerop e)
57 (format nil "~a" c))
58 ((= e 1)
59 (format nil "~a~a" c v))
60 (t
61 (format nil "~a~a^~a" c v e)))))
The exponent type specifier (integer 0 *) in the struct definition is not merely documentation: Common Lisp compilers can use it to generate better code, and the runtime will signal a type error if anything other than a non-negative integer is stored there. The term-negate and term-scale functions are pure: they build new sym-term structs rather than mutating the original. This purity matters because the polynomial arithmetic functions rely on these helpers and expect their inputs to remain unchanged.
term->string uses a three-branch cond to handle the three visually distinct cases: a constant term (exponent zero, print only the number), a linear term (exponent one, omit the ^1), and a higher-degree term (print the full c v^e form). The symbol-name call converts the Lisp symbol stored in the variable to a plain string so that format does not print an unexpected package prefix.
Polynomials
A polynomial is an ordered collection of monomials in a single variable. The implementation stores terms in descending exponent order and automatically combines like terms on every construction. This canonical form means that two polynomials representing the same mathematical object will always have the same in-memory representation.
1 (defstruct (sym-polynomial (:conc-name polynomial-))
2 "A polynomial in one VARIABLE represented as an ordered list of TERMS.
3 Terms are kept in descending order by exponent.
4 DOMAIN is :real (default) or :complex."
5 (variable nil) ; sym-variable
6 (terms '() :type list) ; list of sym-term, descending exponent
7 (domain :real :type keyword))
8
9 (defun %sort-and-combine-terms (terms)
10 "Internal: combine like-exponent terms, drop zero coefficients, sort descending."
11 (let ((table (make-hash-table :test #'eql)))
12 ;; accumulate coefficients by exponent
13 (dolist (term terms)
14 (let ((e (term-exponent term)))
15 (setf (gethash e table)
16 (+ (gethash e table 0)
17 (term-coefficient term)))))
18 ;; rebuild term list, drop zeros, sort descending
19 (let ((result '()))
20 (maphash (lambda (exp coeff)
21 (unless (zerop coeff)
22 (push (cons exp coeff) result)))
23 table)
24 result)))
25
26 (defun make-polynomial (variable terms &key (domain :real))
27 "Create a polynomial in VARIABLE from a list of TERMS.
28 Like-exponent terms are combined automatically; zero-coefficient terms are dropped.
29 Result terms are stored in descending exponent order."
30 (check-type variable sym-variable)
31 (dolist (term terms)
32 (check-type term sym-term)
33 (assert (variable= variable (term-variable term))
34 (term)
35 "Term variable ~s does not match polynomial variable ~s"
36 (variable-name (term-variable term))
37 (variable-name variable)))
38 (assert (member domain '(:real :complex))
39 (domain) "DOMAIN must be :real or :complex, got ~s" domain)
40 (let* ((pairs (%sort-and-combine-terms terms))
41 (sorted (sort pairs #'> :key #'car))
42 (new-terms (mapcar (lambda (pair)
43 (make-sym-term :coefficient (cdr pair)
44 :variable variable
45 :exponent (car pair)))
46 sorted)))
47 (make-sym-polynomial :variable variable
48 :terms new-terms
49 :domain domain)))
50
51 (defun polynomial-degree (poly)
52 "Return the degree (highest exponent) of POLY.
53 Returns -1 for the zero polynomial (no terms)."
54 (check-type poly sym-polynomial)
55 (if (null (polynomial-terms poly))
56 -1
57 (term-exponent (first (polynomial-terms poly)))))
58
59 (defun polynomial-add (p q)
60 "Return P + Q as a new polynomial."
61 (check-type p sym-polynomial)
62 (check-type q sym-polynomial)
63 (assert (variable= (polynomial-variable p) (polynomial-variable q))
64 () "Cannot add polynomials in different variables: ~s and ~s"
65 (variable-name (polynomial-variable p))
66 (variable-name (polynomial-variable q)))
67 (make-polynomial (polynomial-variable p)
68 (append (polynomial-terms p) (polynomial-terms q))
69 :domain (polynomial-domain p)))
70
71 (defun polynomial-negate (poly)
72 "Return -POLY as a new polynomial."
73 (check-type poly sym-polynomial)
74 (make-polynomial (polynomial-variable poly)
75 (mapcar #'term-negate (polynomial-terms poly))
76 :domain (polynomial-domain poly)))
77
78 (defun polynomial-subtract (p q)
79 "Return P - Q as a new polynomial."
80 (polynomial-add p (polynomial-negate q)))
81
82 (defun polynomial-scale (poly scalar)
83 "Return SCALAR * POLY as a new polynomial."
84 (check-type poly sym-polynomial)
85 (check-type scalar real)
86 (make-polynomial (polynomial-variable poly)
87 (mapcar (lambda (term) (term-scale term scalar))
88 (polynomial-terms poly))
89 :domain (polynomial-domain poly)))
90
91 (defun polynomial-evaluate (poly value)
92 "Evaluate POLY at VALUE by substituting the variable."
93 (check-type poly sym-polynomial)
94 (check-type value real)
95 (reduce #'+
96 (polynomial-terms poly)
97 :key (lambda (term)
98 (* (term-coefficient term)
99 (expt value (term-exponent term))))
100 :initial-value 0))
101
102 (defun polynomial->string (poly)
103 "Return a human-readable string for POLY, e.g. '3x^2 + -1x^1 + 5'."
104 (check-type poly sym-polynomial)
105 (if (null (polynomial-terms poly))
106 "0"
107 (format nil "~{~a~^ + ~}"
108 (mapcar #'term->string (polynomial-terms poly)))))
The private helper %sort-and-combine-terms performs the normalization work in two passes. The first pass accumulates coefficients into a hash table keyed by exponent, so any duplicate exponents are automatically merged. The second pass walks the table, drops entries whose total coefficient is zero, and collects the survivors as (exponent . coefficient) cons cells. The public make-polynomial then sorts those pairs in descending order and rebuilds proper sym-term structs.
Polynomial arithmetic is elegantly simple because make-polynomial already normalizes. Addition merely appends the two term lists and lets the constructor handle combination and sorting. Subtraction negates q (by mapping term-negate over its terms) then calls addition. Scaling maps term-scale over all terms. Evaluation uses reduce with an accumulating + over all terms, computing each term’s contribution as c * value^n with Common Lisp’s built-in expt.
Integrals as Data
The last data structure in the layer represents an unevaluated integral expression, either indefinite (no bounds) or definite (with a lower and upper bound). Storing the integral as data rather than immediately computing it allows the system to display it symbolically before choosing to evaluate it numerically.
1 (defstruct (sym-integral (:conc-name integral-))
2 "Represents a definite or indefinite integral.
3 INTEGRAND — the expression being integrated
4 VARIABLE — the variable of integration
5 LOWER — lower bound (number, sym-constant, or NIL for indefinite)
6 UPPER — upper bound (number, sym-constant, or NIL for indefinite)"
7 (integrand nil)
8 (variable nil)
9 (lower nil)
10 (upper nil))
11
12 (defun make-integral (integrand variable &key lower upper)
13 "Create a symbolic integral.
14 Either both LOWER and UPPER must be provided, or neither."
15 (check-type variable sym-variable)
16 (assert (or (and (null lower) (null upper))
17 (and lower upper))
18 () "Either both LOWER and UPPER must be provided, or neither.")
19 (make-sym-integral :integrand integrand
20 :variable variable
21 :lower lower
22 :upper upper))
23
24 (defun integral-definite-p (integral)
25 "Return T if INTEGRAL is a definite integral (has bounds)."
26 (check-type integral sym-integral)
27 (not (null (integral-lower integral))))
28
29 (defun %bound->string (bound)
30 "Internal: render a bound (number or sym-constant) as a string."
31 (cond
32 ((null bound) "")
33 ((numberp bound) (format nil "~a" bound))
34 ((sym-constant-p bound) (symbol-name (constant-name bound)))
35 (t (format nil "~a" bound))))
36
37 (defun integral->string (integral)
38 "Return a human-readable Unicode string for INTEGRAL.
39 Indefinite: ∫(expr) dx
40 Definite: ∫[a,b](expr) dx"
41 (check-type integral sym-integral)
42 (let* ((var (symbol-name (variable-name (integral-variable integral))))
43 (body (cond
44 ((sym-polynomial-p (integral-integrand integral))
45 (polynomial->string (integral-integrand integral)))
46 (t
47 (format nil "~a" (integral-integrand integral)))))
48 (bounds (if (integral-definite-p integral)
49 (format nil "[~a,~a]"
50 (%bound->string (integral-lower integral))
51 (%bound->string (integral-upper integral)))
52 "")))
53 (format nil "∫~a(~a) d~a" bounds body var)))
The make-integral constructor enforces the constraint that bounds must come in pairs: you may supply both lower and upper, or neither, but supplying only one is a programming error caught at runtime with assert. This is a small but valuable design choice because it prevents silently constructing a malformed integral whose lower bound is non-nil but whose upper bound is nil.
integral->string uses the Unicode integral sign ∫ directly in a format string, which works in any modern Common Lisp environment whose source files are saved in UTF-8. The body of the integral is rendered by dispatching on the integrand’s type: if it is a sym-polynomial the existing polynomial->string is used; anything else falls back to format’s default ~a printer. Bounds are rendered by the private %bound->string helper, which knows how to convert both plain numbers and sym-constant objects to strings.
Data Layer Smoke Test
The file closes with a self-contained smoke test that exercises all five types together and prints human-readable results to standard output.
1 (defun run-smoke-test ()
2 "Run a quick sanity check printing results to *standard-output*."
3 (let* ((x (make-variable 'x))
4 (pi-const (make-constant 'pi :pi))
5 ;; Build 3x² - x + 5
6 (p (make-polynomial x
7 (list (make-term 3 x 2)
8 (make-term -1 x 1)
9 (make-term 5 x 0))))
10 ;; Build x + 2
11 (q (make-polynomial x
12 (list (make-term 1 x 1)
13 (make-term 2 x 0))))
14 (sum (polynomial-add p q))
15 (diff (polynomial-subtract p q))
16 (scaled (polynomial-scale p 2))
17 (indef (make-integral p x))
18 (def (make-integral p x :lower 0 :upper 1))
19 (pi-int (make-integral p x :lower 0 :upper pi-const)))
20
21 (format t "~%=== Symbolic Math Data Layer Smoke Test ===~%~%")
22 (format t "p : ~a~%" (polynomial->string p))
23 (format t "q : ~a~%" (polynomial->string q))
24 (format t "p + q : ~a~%" (polynomial->string sum))
25 (format t "p - q : ~a~%" (polynomial->string diff))
26 (format t "2 * p : ~a~%" (polynomial->string scaled))
27 (format t "degree(p) : ~a~%" (polynomial-degree p))
28 (format t "p(0) : ~a~%" (polynomial-evaluate p 0))
29 (format t "p(1) : ~a~%" (polynomial-evaluate p 1))
30 (format t "p(2) : ~a~%" (polynomial-evaluate p 2))
31 (format t "indef : ~a~%" (integral->string indef))
32 (format t "def : ~a~%" (integral->string def))
33 (format t "pi-int : ~a~%" (integral->string pi-int))
34 (format t "~%===========================================~%")))
The smoke test is written as a single large let* binding so that every intermediate result is named and visible. All variables, polynomials, sums, differences, scaled copies, and integral objects are constructed in the binding clauses; the body is purely a sequence of format calls that print them. This layout makes it trivial to inspect the test in a REPL: evaluate any sub-expression from the binding list and examine the result interactively.
Running (run-smoke-test) after loading data.lisp produces output like the following, confirming that polynomial arithmetic, degree computation, evaluation, and integral display all work correctly:
1 === Symbolic Math Data Layer Smoke Test ===
2
3 p : 3x^2 + -1x^1 + 5
4 q : 1x^1 + 2
5 p + q : 3x^2 + 7
6 p - q : 3x^2 + -2x^1 + 3
7 2 * p : 6x^2 + -2x^1 + 10
8 degree(p) : 2
9 p(0) : 5
10 p(1) : 7
11 p(2) : 15
12 indef : ∫(3x^2 + -1x^1 + 5) dx
13 def : ∫[0,1](3x^2 + -1x^1 + 5) dx
14 pi-int : ∫[0,pi](3x^2 + -1x^1 + 5) dx
15
16 ===========================================
Symbolic Differentiation
With the data layer in place we can build the differentiation engine. Polynomial differentiation is driven by two classical rules: the power rule for individual terms, and the sum rule for multi-term polynomials. Because our polynomial representation is already a list of terms, the sum rule requires no special code and it is a natural consequence of mapping the power rule over every element of the list.
The differentiation package lives in differentiation.lisp and declares itself as SYMBOLIC-MATH/DIFF, using (:use #:cl #:symbolic-math) to inherit the full data layer.
The Power Rule on a Single Term
The power rule states that d/dx(c · xⁿ) = n·c · x^(n−1) for n ≥ 1, and that the derivative of a constant is zero. The differentiate-term function encodes this directly.
1 (defun differentiate-term (term)
2 "Apply the power rule to a single monomial TERM.
3
4 d/dx(c · x^n) = n·c · x^(n-1) for n ≥ 1
5 d/dx(c) = 0 for n = 0 (returns NIL)
6
7 Returns a new sym-term, or NIL when the term differentiates to zero."
8 (check-type term sym-term)
9 (let ((c (term-coefficient term))
10 (v (term-variable term))
11 (n (term-exponent term)))
12 (if (zerop n)
13 nil ; constant term → zero, drop it
14 (make-term (* n c) v (1- n)))))
The implementation is a two-branch if. When the exponent is zero the term is a constant and its derivative is zero; returning nil rather than a zero-coefficient term lets the caller use remove nil to discard it cleanly without an extra zero in the term list. For any positive exponent the function multiplies the coefficient by the exponent (n·c) and decrements the exponent by one (n−1), producing a new sym-term via make-term.
Because the function is pure, it reads three slots from an existing term and constructs a brand-new one and it is trivially testable in isolation. It also composes naturally with mapcar since it is a function from one term to one term-or-nil, which is exactly what polynomial-level differentiation needs.
Differentiating a Polynomial
The differentiate function applies differentiate-term to every term in a polynomial, discards the nil results, and re-assembles the survivors into a new polynomial.
1 (defun differentiate (poly)
2 "Differentiate polynomial POLY with respect to its variable.
3
4 Uses the power rule on each term and the sum rule across terms.
5 Returns a new sym-polynomial (the zero polynomial for a constant input).
6
7 Examples:
8 (differentiate p) ; p = 3x² - x + 5 → 6x - 1
9 (differentiate (zero-polynomial x)) → 0 (zero polynomial)"
10 (check-type poly sym-polynomial)
11 (let* ((var (polynomial-variable poly))
12 (diff-terms (remove nil
13 (mapcar #'differentiate-term
14 (polynomial-terms poly)))))
15 (if (null diff-terms)
16 (zero-polynomial var)
17 (make-polynomial var diff-terms :domain (polynomial-domain poly)))))
The sum rule is implicit: mapcar applies differentiate-term independently to every term, and remove nil strips out the constant-term derivatives that returned nil. If all terms were constant the surviving list is empty and zero-polynomial is returned explicitly rather than passing an empty list to make-polynomial. This edge case matters: differentiating a constant polynomial must yield the zero polynomial, and spelling it out explicitly is clearer than relying on make-polynomial’s behaviour for an empty input.
Higher-Order Derivatives
Taking repeated derivatives is expressed as a simple recursion. differentiate-n applies differentiate n times by peeling one application off per recursive call.
1 (defun differentiate-n (poly n)
2 "Return the N-th derivative of polynomial POLY.
3 N must be a non-negative integer.
4
5 (differentiate-n p 0) → p (identity)
6 (differentiate-n p 1) → p'
7 (differentiate-n p 2) → p''
8 …"
9 (check-type poly sym-polynomial)
10 (check-type n (integer 0 *))
11 (cond
12 ((zerop n) poly)
13 (t (differentiate-n (differentiate poly) (1- n)))))
The base case (zerop n) returns the polynomial unchanged, encoding the mathematical identity that the zeroth derivative of a function is the function itself. The recursive case differentiates once and recurses with n-1. Because polynomials have finite degree, the recursion always terminates: at most (degree poly + 1) calls are needed before the result is the zero polynomial, and further applications of differentiate to the zero polynomial keep returning the zero polynomial.
Numerical Gradient and Critical Points
Two convenience functions bridge the symbolic and numeric worlds. gradient-at differentiates a polynomial symbolically and then evaluates the result at a specific point, yielding the slope of the curve at that x-value. critical-point-p uses gradient-at to test whether the derivative is numerically zero at a given point.
1 (defun gradient-at (poly value)
2 "Return the numerical value of the derivative of POLY at VALUE.
3 Differentiates symbolically then evaluates at VALUE.
4
5 (gradient-at p 1) ; slope of p at x=1"
6 (check-type poly sym-polynomial)
7 (check-type value real)
8 (polynomial-evaluate (differentiate poly) value))
9
10 (defun critical-point-p (poly value &key (tolerance 1.0d-9))
11 "Return T if VALUE is (numerically) a critical point of POLY,
12 i.e. |p'(VALUE)| < TOLERANCE.
13
14 (critical-point-p p 1/6) ; is x=1/6 a critical point of p?"
15 (< (abs (gradient-at poly value)) tolerance))
gradient-at is a simple composition: differentiate then evaluate. The symbolic step ensures the derivative is exact; the numeric step converts the exact result to a number. critical-point-p compares the absolute value of the gradient against a configurable tolerance, defaulting to 1.0d-9. Using a tolerance rather than testing for exact equality is important because floating-point arithmetic may introduce small rounding errors even when the true derivative is zero. The &key parameter makes the tolerance easy to adjust without changing the call site.
Differentiation Smoke Test
The differentiation module includes its own smoke test to verify all four exported functions.
1 (defun run-smoke-test ()
2 "Run a quick sanity check on the differentiation functions, printing results."
3 (let* ((x (make-variable 'x))
4
5 ;; p = 3x² - x + 5
6 (p (make-polynomial x (list (make-term 3 x 2)
7 (make-term -1 x 1)
8 (make-term 5 x 0))))
9 (dp (differentiate p)) ; p' = 6x - 1
10 (ddp (differentiate dp)) ; p'' = 6
11 (dddp (differentiate ddp)) ; p''' = 0
12
13 ;; q = x^4 - 2x^3 + x
14 (q (make-polynomial x (list (make-term 1 x 4)
15 (make-term -2 x 3)
16 (make-term 1 x 1))))
17 (dq (differentiate q)) ; q' = 4x^3 - 6x^2 + 1
18 (ddq (differentiate-n q 2)) ; q'' = 12x^2 - 12x
19
20 ;; Critical point of p: p'(x)=0 → 6x-1=0 → x=1/6
21 (cp (/ 1 6)))
22
23 (format t "~%=== Differentiation Smoke Test ===~%~%")
24 (format t "p : ~a~%" (polynomial->string p))
25 (format t "p' : ~a~%" (polynomial->string dp))
26 (format t "p'' : ~a~%" (polynomial->string ddp))
27 (format t "p''' : ~a~%" (polynomial->string dddp))
28 (format t "~%")
29 (format t "q : ~a~%" (polynomial->string q))
30 (format t "q' : ~a~%" (polynomial->string dq))
31 (format t "q'' (via n=2): ~a~%" (polynomial->string ddq))
32 (format t "~%")
33 (format t "gradient-at(p, 0) : ~a (expected -1)~%"
34 (gradient-at p 0))
35 (format t "gradient-at(p, 1) : ~a (expected 5)~%"
36 (gradient-at p 1))
37 (format t "critical-point-p(p,1/6): ~a (expected T)~%"
38 (critical-point-p p cp))
39 (format t "critical-point-p(p,0) : ~a (expected NIL)~%"
40 (critical-point-p p 0))
41 (format t "~%==================================~%")))
The test covers differentiation of a quadratic through to the zero polynomial, a degree-4 polynomial differentiated both step-by-step and with differentiate-n, and two critical-point tests, one at the true critical point x = 1/6 and one at x = 0 which is not a critical point. Embedding the expected values in the format strings (expected -1, expected T, etc.) means the output is self-documenting: a developer reading the terminal output can immediately see whether each result is correct without consulting a separate specification.
Here we run the smoke test:
1 $ sbcl
2 * (load "differentiation.lisp")
3 * (in-package #:symbolic-math/diff)
4 #<package "SYMBOLIC-MATH/DIFF">
5 * (run-smoke-test)
6
7 === Differentiation Smoke Test ===
8
9 p : 3X^2 + -1X + 5
10 p' : 6X + -1
11 p'' : 6
12 p''' : 0
13
14 q : 1X^4 + -2X^3 + 1X
15 q' : 4X^3 + -6X^2 + 1
16 q'' (via n=2): 12X^2 + -12X
17
18 gradient-at(p, 0) : -1 (expected -1)
19 gradient-at(p, 1) : 5 (expected 5)
20 critical-point-p(p,1/6): t (expected T)
21 critical-point-p(p,0) : nil (expected NIL)
22
23 ==================================
24 nil
25 *
Symbolic Integration
Integration is the inverse of differentiation. Where differentiation reduces degree by one, integration raises it. The integration package in integration.lisp implements the reverse power rule for individual terms, extends it to polynomials via the sum rule, evaluates definite integrals numerically using the Fundamental Theorem of Calculus, and provides helpers for constructing sym-integral shell objects.
The package is SYMBOLIC-MATH/INTEG, again using (:use #:cl #:symbolic-math).
The Reverse Power Rule on a Single Term
The reverse power rule states that ∫ c · xⁿ dx = (c / (n+1)) · x^(n+1). Because n is always a non-negative integer here, n+1 is always positive, so division is always well-defined and no special case is needed.
1 (defun integrate-term (term)
2 "Apply the reverse power rule to a single monomial TERM.
3
4 ∫ c · x^n dx = (c / (n+1)) · x^(n+1)
5
6 Returns a new sym-term. The constant of integration (+C) is handled at
7 the polynomial level by the caller.
8
9 Examples:
10 (integrate-term (make-term 3 x 2)) ; => (1)x^3 [3/(2+1) = 1]
11 (integrate-term (make-term -1 x 1)) ; => (-1/2)x^2
12 (integrate-term (make-term 5 x 0)) ; => 5x^1"
13 (check-type term sym-term)
14 (let* ((c (term-coefficient term))
15 (v (term-variable term))
16 (n (term-exponent term))
17 (new-exp (1+ n))
18 (new-coef (/ c new-exp))) ; exact ratio; CL keeps it rational
19 (make-term new-coef v new-exp)))
The division (/ c new-exp) is the algebraic heart of the function. In Common Lisp, dividing an integer by another integer produces an exact rational number rather than a float. Integrating the term (make-term 3 x 2) (representing 3x²) yields the new coefficient 3/3 = 1; integrating (make-term -1 x 1) (representing −x) yields -1/2. These rationals propagate unchanged through subsequent operations, so no precision is lost until the caller explicitly requests a float which is a significant advantage over languages that default to floating-point division.
Unlike differentiate-term, this function never returns nil because every term, including a constant term with exponent zero, integrates to a non-zero result. A constant 5 integrates to 5x, so there is no zero-result case to handle.
The Antiderivative of a Polynomial
The function integrate maps integrate-term over every term of the polynomial and assembles the results into a new polynomial of degree (degree poly + 1).
1 (defun integrate (poly)
2 "Return the antiderivative of polynomial POLY (without the constant of
3 integration +C). Applies the reverse power rule to each term.
4
5 Examples:
6 poly = 3x² - x + 5
7 ∫poly = x³ - (1/2)x² + 5x
8
9 Returns a new sym-polynomial of degree (degree(POLY) + 1)."
10 (check-type poly sym-polynomial)
11 (let* ((var (polynomial-variable poly))
12 (new-terms (mapcar #'integrate-term (polynomial-terms poly))))
13 (if (null new-terms)
14 (zero-polynomial var)
15 (make-polynomial var new-terms :domain (polynomial-domain poly)))))
The structure mirrors differentiate almost exactly, with integrate-term replacing differentiate-term and remove nil omitted because integrate-term never returns nil. The zero-polynomial edge case is still handled: integrating the zero polynomial (no terms) returns the zero polynomial rather than passing an empty list to make-polynomial. Notice that the constant of integration +C is deliberately omitted. In an indefinite integration context +C represents an entire family of functions; adding a zero-exponent term for it would require an arbitrary choice for the constant’s coefficient. By leaving it out, the function returns the canonical antiderivative, and callers can add their own constant term if needed.
Definite Integral Evaluation
The Fundamental Theorem of Calculus reduces a definite integral to two polynomial evaluations: ∫[a,b] f dx = F(b) − F(a), where F is the antiderivative of f.
1 (defun %resolve-bound (bound)
2 "Internal: convert a bound (number or sym-constant) to a number."
3 (cond
4 ((null bound)
5 (error "Bound is NIL — use evaluate-definite only on definite integrals."))
6 ((numberp bound) bound)
7 ((constant-p bound) (constant-numeric-value bound))
8 (t (error "Unrecognised bound type: ~s" bound))))
9
10 (defun evaluate-definite (poly lower upper)
11 "Evaluate the definite integral ∫[LOWER,UPPER] POLY dx numerically.
12
13 LOWER and UPPER may be real numbers or sym-constant objects.
14
15 Uses the Fundamental Theorem of Calculus:
16 ∫[a,b] f dx = F(b) - F(a) where F = (integrate poly)
17
18 Examples:
19 (evaluate-definite p 0 1)
20 (evaluate-definite p 0 (make-constant 'pi :pi))"
21 (check-type poly sym-polynomial)
22 (let* ((F (integrate poly))
23 (a (coerce (%resolve-bound lower) 'double-float))
24 (b (coerce (%resolve-bound upper) 'double-float)))
25 (- (polynomial-evaluate F b)
26 (polynomial-evaluate F a))))
The function %resolve-bound is a small dispatch function that converts whichever bound type was supplied: a plain number, sym-constant, or an unexpected type to a numeric value or signals an informative error. The two coerce calls in evaluate-definite ensure the bounds are double-float before passing them to polynomial-evaluate, which prevents mixed exact/floating-point arithmetic from producing surprising results. The final subtraction F(b) - F(a) is the Fundamental Theorem in code, and its simplicity is a direct reflection of the theorem’s elegance.
Iterated Antiderivatives
Just as differentiate-n applies differentiation n times, integrate-n applies integration n times.
1 (defun integrate-n (poly n)
2 "Return the N-th iterated antiderivative of POLY.
3 N must be a non-negative integer.
4
5 (integrate-n p 0) → p (identity)
6 (integrate-n p 1) → ∫p dx
7 (integrate-n p 2) → ∫∫p dx dx
8 …"
9 (check-type poly sym-polynomial)
10 (check-type n (integer 0 *))
11 (cond
12 ((zerop n) poly)
13 (t (integrate-n (integrate poly) (1- n)))))
The recursion mirrors differentiate-n in structure: base case returns the polynomial unchanged; recursive case integrates once and decrements n. Because each application of integrate raises the degree by exactly one, integrate-n on a polynomial of degree d produces a polynomial of degree d + n. The rational arithmetic in integrate-term becomes especially valuable here: integrating twice starting from integer coefficients yields rational coefficients that are still exact, whereas floating-point arithmetic would accumulate rounding error with each pass.
Integral Shell Constructors
Two thin wrappers make it convenient to package a polynomial together with its bounds into the sym-integral data structure defined in the data layer.
1 (defun make-indefinite-integral (poly)
2 "Create a sym-integral wrapping POLY (indefinite form).
3 The stored integrand is POLY itself.
4
5 Rendering: ∫(expr) dx"
6 (check-type poly sym-polynomial)
7 (make-integral poly (polynomial-variable poly)))
8
9 (defun make-definite-integral (poly lower upper)
10 "Create a sym-integral wrapping POLY with explicit bounds LOWER and UPPER.
11 LOWER and UPPER may be real numbers or sym-constant objects.
12
13 Rendering: ∫[a,b](expr) dx"
14 (check-type poly sym-polynomial)
15 (make-integral poly (polynomial-variable poly) :lower lower :upper upper))
These functions exist purely for ergonomics. Without them, a caller would need to separately extract the polynomial’s variable and pass it as the second argument to make-integral that is a minor but unnecessary repetition. By encapsulating that extraction, the wrappers let calling code read more naturally: (make-indefinite-integral p) rather than (make-integral p (polynomial-variable p)). Note that these constructors store the original polynomial as the integrand; they do not pre-compute the antiderivative. Evaluation is deferred to an explicit call to evaluate-definite or integrate.
Integration Smoke Test
The integration smoke test verifies each exported function and echoes expected values inline.
1 (defun run-smoke-test ()
2 "Run a quick sanity check on the integration functions, printing results."
3 (let* ((x (make-variable 'x))
4 (pi-c (make-constant 'pi :pi))
5
6 ;; p = 3x² - x + 5
7 (p (make-polynomial x (list (make-term 3 x 2)
8 (make-term -1 x 1)
9 (make-term 5 x 0))))
10 (ip (integrate p)) ; x³ - (1/2)x² + 5x
11 (iip (integrate-n p 2)) ; (1/4)x⁴ - (1/6)x³ + (5/2)x²
12
13 ;; q = 6x + 2
14 (q (make-polynomial x (list (make-term 6 x 1)
15 (make-term 2 x 0))))
16 (iq (integrate q)) ; 3x² + 2x
17
18 (def-01 (evaluate-definite p 0 1)) ; expected 5.5
19 (def-q-01 (evaluate-definite q 0 1)) ; expected 5.0
20 (def-pi (evaluate-definite p 0 pi-c))
21
22 (indef-shell (make-indefinite-integral p))
23 (def-shell (make-definite-integral p 0 1))
24 (pi-shell (make-definite-integral p 0 pi-c)))
25
26 (format t "~%=== Integration Smoke Test ===~%~%")
27 (format t "p : ~a~%" (polynomial->string p))
28 (format t "∫p dx : ~a~%" (polynomial->string ip))
29 (format t "∫∫p dx dx : ~a~%" (polynomial->string iip))
30 (format t "~%")
31 (format t "q : ~a~%" (polynomial->string q))
32 (format t "∫q dx : ~a~%" (polynomial->string iq))
33 (format t "~%")
34 (format t "∫₀¹ p dx : ~a (expected 5.5)~%" def-01)
35 (format t "∫₀¹ q dx : ~a (expected 5.0)~%" def-q-01)
36 (format t "∫₀^π p dx : ~f~%" def-pi)
37 (format t "~%")
38 (format t "indef shell : ~a~%" (integral->string indef-shell))
39 (format t "def [0,1] : ~a~%" (integral->string def-shell))
40 (format t "def [0,pi] : ~a~%" (integral->string pi-shell))
41 (format t "~%==============================~%")))
The test checks the antiderivative of p (a quadratic), the second antiderivative of p, the antiderivative of q (a linear polynomial), two numerical definite integrals with known exact values (5.5 and 5.0), a definite integral with a sym-constant bound (π), and the string rendering of all three integral shell variants. The ~f directive for the π-bounded result prints a floating-point number whose exact digits will vary by implementation but whose value should be approximately 102.olean that is a sanity check that the numeric path through sym-constant works end-to-end.
Here is the output of the smoke test:
1 $ sbcl
2 * (load "integration.lisp")
3 warning: redefining symbolic-math:run-smoke-test in DEFUN
4 t
5 * (in-package #:symbolic-math/integ)
6 #<package "SYMBOLIC-MATH/INTEG">
7 * (run-smoke-test)
8
9 === Integration Smoke Test ===
10
11 p : 3X^2 + -1X + 5
12 ∫p dx : 1X^3 + -1/2X^2 + 5X
13 ∫∫p dx dx : 1/4X^4 + -1/6X^3 + 5/2X^2
14
15 q : 6X + 2
16 ∫q dx : 3X^2 + 2X
17
18 ∫₀¹ p dx : 5.5d0 (expected 5.5)
19 ∫₀¹ q dx : 5.0d0 (expected 5.0)
20 ∫₀^π p dx : 41.77943774770411
21
22 indef shell : ∫(3X^2 + -1X + 5) dX
23 def [0,1] : ∫[0,1](3X^2 + -1X + 5) dX
24 def [0,pi] : ∫[0,PI](3X^2 + -1X + 5) dX
25
26 ==============================
27 nil
28 *
Wrap Up
In this chapter we built a three-file symbolic mathematics library in Common Lisp that demonstrates how a language designed for symbolic computation can express mathematical concepts cleanly and correctly.
The data layer (data.lisp) showed how defstruct provides a lightweight, functional alternative to CLOS for domain objects. By enforcing invariants at construction time with check-type and assert, we ensured that ill-formed objects never enter the system. The canonical-form polynomial representation: terms always sorted by descending exponent with like terms combined, removing the need for normalization logic in every downstream consumer.
The differentiation layer (differentiation.lisp) demonstrated that the power rule and sum rule can be encoded almost literally in code. differentiate-term is a direct transliteration of d/dx(cxⁿ) = ncx^(n−1); differentiate applies it via mapcar and remove nil, and the sum rule emerges automatically from the polynomial’s list structure. Higher-order derivatives and numerical gradient evaluation required only a handful of additional lines.
The integration layer (integration.lisp) mirrored the differentiation layer but introduced the important advantage of Common Lisp’s exact rational arithmetic. Coefficients produced by the reverse power rule, such as the -1/2 arising from integrating -x and are stored as exact rationals, so no precision is lost across multiple integrations. The Fundamental Theorem of Calculus translated into a three-line function that computes an antiderivative symbolically and then performs two numerical evaluations, achieving both symbolic clarity and numeric accuracy.
Taken together, the three files illustrate a broader lesson: when you model a problem domain faithfully as data, the algorithms that manipulate that data often become nearly self-evident. The mathematics drives the code structure rather than the other way around, and the result is software that is simultaneously easier to understand, easier to test, and easier to extend.