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 following diagram shows the high-level architecture of the system:
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).
Packages
The three files define three packages. Only the data layer uses #:cl alone; the differentiation and integration packages inherit the whole data layer through :use, so they can refer to make-polynomial, polynomial->string, and friends without any package prefix.
1 (defpackage #:symbolic-math
2 (:use #:cl)
3 (:export
4 ;; Variables
5 #:sym-variable
6 #:make-variable
7 #:variable-p
8 #:variable-name
9 #:variable-domain
10 #:variable=
11 ;; Constants
12 #:sym-constant
13 #:make-constant
14 #:constant-p
15 #:constant-name
16 #:constant-value
17 #:constant-numeric-value
18 ;; Terms (monomials)
19 #:sym-term
20 #:make-term
21 #:term-p
22 #:term-coefficient
23 #:term-variable
24 #:term-exponent
25 #:term=
26 #:term-negate
27 #:term-scale
28 #:term->string
29 ;; Polynomials
30 #:sym-polynomial
31 #:make-polynomial
32 #:polynomial-p
33 #:polynomial-variable
34 #:polynomial-terms
35 #:polynomial-domain
36 #:polynomial-degree
37 #:polynomial-leading-term
38 #:polynomial-add
39 #:polynomial-subtract
40 #:polynomial-negate
41 #:polynomial-scale
42 #:polynomial->string
43 #:polynomial-evaluate
44 #:polynomial-normalize
45 ;; Polynomial convenience constructors
46 #:zero-polynomial
47 #:constant-polynomial
48 #:identity-polynomial
49 ;; Integrals
50 #:sym-integral
51 #:make-integral
52 #:integral-p
53 #:integral-integrand
54 #:integral-variable
55 #:integral-lower
56 #:integral-upper
57 #:integral-definite-p
58 #:integral->string
59 ;; Smoke test
60 #:run-smoke-test))
Both of the dependent packages shadow run-smoke-test, and the reason is worth pausing over because it is a classic Common Lisp trap. Because SYMBOLIC-MATH exports run-smoke-test and the other two packages :use it, that symbol is accessible in SYMBOLIC-MATH/DIFF and SYMBOLIC-MATH/INTEG without being present in them. A bare (defun run-smoke-test ...) would therefore resolve to the inherited SYMBOLIC-MATH:RUN-SMOKE-TEST and redefine the data layer’s smoke test rather than define a new one. All three modules would then share a single function, and the last file loaded would win. Adding (:shadow #:run-smoke-test) to the package definition forces a fresh, package-local symbol so that each module really does have its own test.
The library also ships an ASDF system definition, symbolic-math.asd, so it can be loaded with (asdf:load-system :symbolic-math) or (ql:quickload :symbolic-math) instead of loading the three files by hand:
1 (asdf:defsystem #:symbolic-math
2 :description "A small, purely functional symbolic mathematics library for Common Lisp: variables, constants, monomials, polynomials, symbolic differentiation, and symbolic integration."
3 :author "Mark Watson"
4 :license "Apache-2.0"
5 :serial t
6 :components ((:file "data")
7 (:file "differentiation")
8 (:file "integration")))
The :serial t guarantees that data.lisp is compiled and loaded before the two layers that depend on it. When you load differentiation.lisp or integration.lisp directly instead, the (unless (find-package '#:symbolic-math) (load ...)) form at the top of each file loads data.lisp first, so both workflows work.
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 (check-type name symbol)
15 (assert (member domain '(:real :complex :integer))
16 (domain) "DOMAIN must be :real, :complex, or :integer, got ~s" domain)
17 (make-sym-variable :name name :domain domain))
18
19 (defun variable-p (obj)
20 "Return T if OBJ is a sym-variable."
21 (sym-variable-p obj))
22
23 (defun variable= (a b)
24 "Return T if variables A and B represent the same symbolic variable.
25 Two variables are equal when their names and domains match."
26 (and (sym-variable-p a)
27 (sym-variable-p b)
28 (eq (variable-name a) (variable-name b))
29 (eq (variable-domain a) (variable-domain b))))
30
31 (defstruct (sym-constant (:conc-name constant-))
32 "A named mathematical constant.
33 VALUE may be a Common Lisp number, :pi, or :e (Euler's number)."
34 (name nil :type symbol)
35 (value nil)) ; number | :pi | :e
36
37 (defun make-constant (name value)
38 "Create a named constant with NAME (symbol) and VALUE.
39 VALUE may be a real number, :pi, or :e.
40
41 Examples:
42 (make-constant 'pi :pi)
43 (make-constant 'e :e)
44 (make-constant 'g 9.80665)"
45 (check-type name symbol)
46 (assert (or (numberp value) (member value '(:pi :e)))
47 (value) "VALUE must be a number, :pi, or :e, got ~s" value)
48 (make-sym-constant :name name :value value))
49
50 (defun constant-p (obj)
51 "Return T if OBJ is a sym-constant."
52 (sym-constant-p obj))
53
54 (defun constant-numeric-value (c)
55 "Return the numeric (floating-point) value of constant C.
56 :pi → pi, :e → e, numbers are returned as-is."
57 (check-type c sym-constant)
58 (let ((v (constant-value c)))
59 (cond
60 ((eq v :pi) pi)
61 ((eq v :e) (exp 1.0d0))
62 (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 the two keyword cases to a double-float on demand using a short cond dispatch, and returns a stored number unchanged. Keeping the symbolic value alongside the name means that bounds such as π can be displayed with their mathematical glyph 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 — 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 (variable ignored for evaluation)"
16 (check-type coefficient real)
17 (check-type variable sym-variable)
18 (check-type exponent (integer 0 *))
19 (make-sym-term :coefficient coefficient
20 :variable variable
21 :exponent exponent))
22
23 (defun term-p (obj)
24 "Return T if OBJ is a sym-term."
25 (sym-term-p obj))
26
27 (defun term= (a b)
28 "Return T if terms A and B are structurally equal
29 (same variable, exponent, and coefficient)."
30 (and (sym-term-p a) (sym-term-p b)
31 (= (term-coefficient a) (term-coefficient b))
32 (variable= (term-variable a) (term-variable b))
33 (= (term-exponent a) (term-exponent b))))
34
35 (defun term-negate (term)
36 "Return a new term equal to -TERM."
37 (make-sym-term :coefficient (- (term-coefficient term))
38 :variable (term-variable term)
39 :exponent (term-exponent term)))
40
41 (defun term-scale (term scalar)
42 "Return a new term equal to SCALAR * TERM."
43 (check-type scalar real)
44 (make-sym-term :coefficient (* scalar (term-coefficient term))
45 :variable (term-variable term)
46 :exponent (term-exponent term)))
47
48 (defun term->string (term)
49 "Return a human-readable string representation of TERM.
50 The variable name is down-cased, so the display matches the mathematical
51 notation used in this library's documentation.
52
53 Example: 3x^2, -1x, 5 (for exponent 0)."
54 (let ((c (term-coefficient term))
55 (v (string-downcase (symbol-name (variable-name (term-variable term)))))
56 (e (term-exponent term)))
57 (cond
58 ((zerop e)
59 (format nil "~a" c))
60 ((= e 1)
61 (format nil "~a~a" c v))
62 (t
63 (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 variable name is printed through (string-downcase (symbol-name ...)), which suppresses any package prefix that format might otherwise emit and also avoids the upper-case X that symbol-name alone would produce for a reader symbol. Note that a coefficient of one is not elided: the linear term of x prints as 1x, not x.
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, and return
11 (exponent . coefficient) pairs sorted by descending exponent."
12 (let ((table (make-hash-table :test #'eql)))
13 ;; accumulate coefficients by exponent
14 (dolist (term terms)
15 (let ((e (term-exponent term)))
16 (setf (gethash e table)
17 (+ (gethash e table 0)
18 (term-coefficient term)))))
19 ;; collect the surviving pairs and sort them by descending exponent
20 (let ((result '()))
21 (maphash (lambda (exp coeff)
22 (unless (zerop coeff)
23 (push (cons exp coeff) result)))
24 table)
25 (sort result #'> :key #'car))))
26
27 (defun make-polynomial (variable terms &key (domain :real))
28 "Create a polynomial in VARIABLE from a list of TERMS (sym-term objects).
29 Like-exponent terms are combined automatically; zero-coefficient terms are dropped.
30 Result terms are stored in descending exponent order.
31
32 Example:
33 (let ((x (make-variable 'x)))
34 (make-polynomial x (list (make-term 3 x 2)
35 (make-term -1 x 1)
36 (make-term 5 x 0))))"
37 (check-type variable sym-variable)
38 (dolist (term terms)
39 (check-type term sym-term)
40 (assert (variable= variable (term-variable term))
41 (term)
42 "Term variable ~s does not match polynomial variable ~s"
43 (variable-name (term-variable term))
44 (variable-name variable)))
45 (assert (member domain '(:real :complex))
46 (domain) "DOMAIN must be :real or :complex, got ~s" domain)
47 (let* ((pairs (%sort-and-combine-terms terms))
48 (new-terms (mapcar (lambda (pair)
49 (make-sym-term :coefficient (cdr pair)
50 :variable variable
51 :exponent (car pair)))
52 pairs)))
53 (make-sym-polynomial :variable variable
54 :terms new-terms
55 :domain domain)))
56
57 (defun polynomial-degree (poly)
58 "Return the degree (highest exponent) of POLY.
59 Returns -1 for the zero polynomial (no terms)."
60 (check-type poly sym-polynomial)
61 (if (null (polynomial-terms poly))
62 -1
63 (term-exponent (first (polynomial-terms poly)))))
64
65 (defun %join-domains (a b)
66 "Internal: return the wider of domains A and B (:real ⊂ :complex).
67 Arithmetic on polynomials from different domains happens in the wider one:
68 a real polynomial plus a complex polynomial is a complex polynomial."
69 (if (or (eq a :complex) (eq b :complex)) :complex :real))
70
71 (defun polynomial-add (p q)
72 "Return P + Q as a new polynomial.
73 P and Q must share the same variable; the result takes the wider of the
74 two domains."
75 (check-type p sym-polynomial)
76 (check-type q sym-polynomial)
77 (assert (variable= (polynomial-variable p) (polynomial-variable q))
78 () "Cannot add polynomials in different variables: ~s and ~s"
79 (variable-name (polynomial-variable p))
80 (variable-name (polynomial-variable q)))
81 (make-polynomial (polynomial-variable p)
82 (append (polynomial-terms p) (polynomial-terms q))
83 :domain (%join-domains (polynomial-domain p)
84 (polynomial-domain q))))
85
86 (defun polynomial-negate (poly)
87 "Return -POLY as a new polynomial."
88 (check-type poly sym-polynomial)
89 (make-polynomial (polynomial-variable poly)
90 (mapcar #'term-negate (polynomial-terms poly))
91 :domain (polynomial-domain poly)))
92
93 (defun polynomial-subtract (p q)
94 "Return P - Q as a new polynomial."
95 (polynomial-add p (polynomial-negate q)))
96
97 (defun polynomial-scale (poly scalar)
98 "Return SCALAR * POLY as a new polynomial."
99 (check-type poly sym-polynomial)
100 (check-type scalar real)
101 (make-polynomial (polynomial-variable poly)
102 (mapcar (lambda (term) (term-scale term scalar))
103 (polynomial-terms poly))
104 :domain (polynomial-domain poly)))
105
106 (defun polynomial-evaluate (poly value)
107 "Evaluate POLY at VALUE (a number) by substituting the variable.
108 Returns a number.
109
110 Example: (polynomial-evaluate p 2) ; evaluate p(2)"
111 (check-type poly sym-polynomial)
112 (check-type value real)
113 (reduce #'+
114 (polynomial-terms poly)
115 :key (lambda (term)
116 (* (term-coefficient term)
117 (expt value (term-exponent term))))
118 :initial-value 0))
119
120 (defun polynomial->string (poly)
121 "Return a human-readable string representation of POLY.
122 Terms are printed in descending exponent order, separated by ' + '.
123 The zero polynomial is rendered as '0'.
124
125 Example: '3x^2 + -1x + 5'"
126 (check-type poly sym-polynomial)
127 (if (null (polynomial-terms poly))
128 "0"
129 (format nil "~{~a~^ + ~}"
130 (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, collects the survivors as (exponent . coefficient) cons cells, and sorts them into descending exponent order. Sorting has to happen here rather than in make-polynomial, because maphash visits entries in an unspecified order: without an explicit sort the pairs would come back out of the table in essentially arbitrary order. The public make-polynomial then rebuilds proper sym-term structs from those already-sorted pairs.
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. It does have to reconcile one piece of metadata, though: the two operands may have different domains, so %join-domains returns the wider of the pair (:real ⊂ :complex) and the sum is built in that domain. A real polynomial plus a complex polynomial is a complex polynomial, and silently inheriting the left operand’s domain would quietly discard that fact. Subtraction negates q (by mapping term-negate over its terms) and then calls addition, so it inherits the same handling. 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
4 Fields:
5 INTEGRAND — the expression being integrated (a sym-polynomial or any form)
6 VARIABLE — the variable of integration (sym-variable)
7 LOWER — lower bound (a number, sym-constant, or NIL for indefinite)
8 UPPER — upper bound (a number, sym-constant, or NIL for indefinite)"
9 (integrand nil) ; sym-polynomial | expression
10 (variable nil) ; sym-variable
11 (lower nil) ; number | sym-constant | nil
12 (upper nil)) ; number | sym-constant | nil
13
14 (defun make-integral (integrand variable &key lower upper)
15 "Create a symbolic integral.
16
17 INTEGRAND — the expression to integrate (typically a sym-polynomial)
18 VARIABLE — sym-variable: the variable of integration
19 LOWER — optional lower bound (number or sym-constant); NIL ⇒ indefinite
20 UPPER — optional upper bound (number or sym-constant); NIL ⇒ indefinite
21
22 Examples:
23 ;; Indefinite: ∫ (3x² - x + 5) dx
24 (make-integral poly x)
25
26 ;; Definite: ∫₀¹ (3x² - x + 5) dx
27 (make-integral poly x :lower 0 :upper 1)"
28 (check-type variable sym-variable)
29 (assert (or (and (null lower) (null upper))
30 (and lower upper))
31 () "Either both LOWER and UPPER must be provided, or neither.")
32 (make-sym-integral :integrand integrand
33 :variable variable
34 :lower lower
35 :upper upper))
36
37 (defun integral-definite-p (integral)
38 "Return T if INTEGRAL is a definite integral (has bounds)."
39 (check-type integral sym-integral)
40 (not (null (integral-lower integral))))
41
42 (defun %bound->string (bound)
43 "Internal: render a bound (number or sym-constant) as a string.
44 Symbolic constants are rendered with their mathematical glyph where one
45 exists (:pi → π, :e → e) and otherwise by their down-cased name."
46 (cond
47 ((null bound) "")
48 ((numberp bound) (format nil "~a" bound))
49 ((sym-constant-p bound)
50 (case (constant-value bound)
51 (:pi "π")
52 (:e "e")
53 (t (string-downcase (symbol-name (constant-name bound))))))
54 (t (format nil "~a" bound))))
55
56 (defun integral->string (integral)
57 "Return a human-readable string for INTEGRAL.
58
59 Indefinite: ∫(expr) dx
60 Definite: ∫[a,b](expr) dx"
61 (check-type integral sym-integral)
62 (let* ((var (string-downcase
63 (symbol-name (variable-name (integral-variable integral)))))
64 (body (cond
65 ((sym-polynomial-p (integral-integrand integral))
66 (polynomial->string (integral-integrand integral)))
67 (t
68 (format nil "~a" (integral-integrand integral)))))
69 (bounds (if (integral-definite-p integral)
70 (format nil "[~a,~a]"
71 (%bound->string (integral-lower integral))
72 (%bound->string (integral-upper integral)))
73 "")))
74 (format nil "∫~a(~a) d~a" bounds body var)))
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 ;; p + q => 3x² + 7
15 (sum (polynomial-add p q))
16 ;; p - q => 3x² - 2x + 3
17 (diff (polynomial-subtract p q))
18 ;; 2p => 6x² - 2x + 10
19 (scaled (polynomial-scale p 2))
20 ;; Indefinite integral of p
21 (indef (make-integral p x))
22 ;; Definite integral ∫₀¹ p dx
23 (def (make-integral p x :lower 0 :upper 1))
24 ;; Definite integral ∫₀^π p dx (π as sym-constant bound)
25 (pi-int (make-integral p x :lower 0 :upper pi-const)))
26
27 (format t "~%=== Symbolic Math Data Layer Smoke Test ===~%~%")
28 (format t "p : ~a~%" (polynomial->string p))
29 (format t "q : ~a~%" (polynomial->string q))
30 (format t "p + q : ~a~%" (polynomial->string sum))
31 (format t "p - q : ~a~%" (polynomial->string diff))
32 (format t "2 * p : ~a~%" (polynomial->string scaled))
33 (format t "degree(p) : ~a~%" (polynomial-degree p))
34 (format t "p(0) : ~a~%" (polynomial-evaluate p 0))
35 (format t "p(1) : ~a~%" (polynomial-evaluate p 1))
36 (format t "p(2) : ~a~%" (polynomial-evaluate p 2))
37 (format t "indef : ~a~%" (integral->string indef))
38 (format t "def : ~a~%" (integral->string def))
39 (format t "pi-int : ~a~%" (integral->string pi-int))
40 (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 $ sbcl
2 * (load "data.lisp")
3 T
4 * (in-package #:symbolic-math)
5 #<PACKAGE "SYMBOLIC-MATH">
6 * (run-smoke-test)
7
8 === Symbolic Math Data Layer Smoke Test ===
9
10 p : 3x^2 + -1x + 5
11 q : 1x + 2
12 p + q : 3x^2 + 7
13 p - q : 3x^2 + -2x + 3
14 2 * p : 6x^2 + -2x + 10
15 degree(p) : 2
16 p(0) : 5
17 p(1) : 7
18 p(2) : 15
19 indef : ∫(3x^2 + -1x + 5) dx
20 def : ∫[0,1](3x^2 + -1x + 5) dx
21 pi-int : ∫[0,π](3x^2 + -1x + 5) dx
22
23 ===========================================
24 NIL
25 *
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:
1 (defpackage #:symbolic-math/diff
2 (:use #:cl #:symbolic-math)
3 ;; SYMBOLIC-MATH exports RUN-SMOKE-TEST, and :USE makes that symbol
4 ;; accessible here. Without this SHADOW, the DEFUN below would redefine
5 ;; SYMBOLIC-MATH:RUN-SMOKE-TEST instead of defining a package-local one.
6 (:shadow #:run-smoke-test)
7 (:export
8 ;; Core differentiation
9 #:differentiate-term
10 #:differentiate
11 ;; Higher-order derivatives
12 #:differentiate-n
13 ;; Convenience
14 #:gradient-at
15 #:critical-point-p
16 ;; Smoke test
17 #:run-smoke-test))
The :shadow line deserves a second look, because without it the module silently breaks. SYMBOLIC-MATH exports run-smoke-test, and :use makes every exported symbol accessible in SYMBOLIC-MATH/DIFF — but accessible is not the same as present. A plain (defun run-smoke-test ...) would resolve to the inherited symbol and redefine SYMBOLIC-MATH:RUN-SMOKE-TEST rather than create a new function. All three files would then define the same symbol, and the last one loaded would win. :shadow interposes a fresh package-local symbol so each file gets its own test.
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, which makes it 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 :domain (polynomial-domain poly))
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. make-polynomial would in fact cope with an empty list, so this branch is about clarity as much as correctness — but it does carry one piece of real information, because zero-polynomial is called with :domain (polynomial-domain poly). Differentiating a complex polynomial must yield a complex zero polynomial, and the default domain of zero-polynomial is :real.
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))
Differentiation Smoke Test
The differentiation module includes its own smoke test to verify all four exported functions. It is genuinely its own: before the :shadow fix described above, this defun redefined the data layer’s smoke test instead of defining a new one.
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
10 ;; p' = 6x - 1
11 (dp (differentiate p))
12
13 ;; p'' = 6
14 (ddp (differentiate dp))
15
16 ;; p''' = 0 (zero polynomial)
17 (dddp (differentiate ddp))
18
19 ;; q = x^4 - 2x^3 + x
20 (q (make-polynomial x (list (make-term 1 x 4)
21 (make-term -2 x 3)
22 (make-term 1 x 1))))
23 ;; q' = 4x^3 - 6x^2 + 1
24 (dq (differentiate q))
25
26 ;; q'' = 12x^2 - 12x
27 (ddq (differentiate-n q 2))
28
29 ;; Critical point of p: p'(x)=0 → 6x-1=0 → x=1/6
30 (cp (/ 1 6)))
31
32 (format t "~%=== Differentiation Smoke Test ===~%~%")
33 (format t "p : ~a~%" (polynomial->string p))
34 (format t "p' : ~a~%" (polynomial->string dp))
35 (format t "p'' : ~a~%" (polynomial->string ddp))
36 (format t "p''' : ~a~%" (polynomial->string dddp))
37 (format t "~%")
38 (format t "q : ~a~%" (polynomial->string q))
39 (format t "q' : ~a~%" (polynomial->string dq))
40 (format t "q'' (via n=2): ~a~%" (polynomial->string ddq))
41 (format t "~%")
42 (format t "gradient-at(p, 0) : ~a (expected -1)~%"
43 (gradient-at p 0))
44 (format t "gradient-at(p, 1) : ~a (expected 5)~%"
45 (gradient-at p 1))
46 (format t "critical-point-p(p,1/6): ~a (expected T)~%"
47 (critical-point-p p cp))
48 (format t "critical-point-p(p,0) : ~a (expected NIL)~%"
49 (critical-point-p p 0))
50 (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 T
4 * (in-package #:symbolic-math/diff)
5 #<PACKAGE "SYMBOLIC-MATH/DIFF">
6 * (run-smoke-test)
7
8 === Differentiation Smoke Test ===
9
10 p : 3x^2 + -1x + 5
11 p' : 6x + -1
12 p'' : 6
13 p''' : 0
14
15 q : 1x^4 + -2x^3 + 1x
16 q' : 4x^3 + -6x^2 + 1
17 q'' (via n=2): 12x^2 + -12x
18
19 gradient-at(p, 0) : -1 (expected -1)
20 gradient-at(p, 1) : 5 (expected 5)
21 critical-point-p(p,1/6): T (expected T)
22 critical-point-p(p,0) : NIL (expected NIL)
23
24 ==================================
25 NIL
26 *
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) and again shadowing run-smoke-test:
1 (defpackage #:symbolic-math/integ
2 (:use #:cl #:symbolic-math)
3 ;; SYMBOLIC-MATH exports RUN-SMOKE-TEST, and :USE makes that symbol
4 ;; accessible here. Without this SHADOW, the DEFUN below would redefine
5 ;; SYMBOLIC-MATH:RUN-SMOKE-TEST instead of defining a package-local one.
6 (:shadow #:run-smoke-test)
7 (:export
8 ;; Core integration
9 #:integrate-term ; reverse power rule on one sym-term
10 #:integrate ; antiderivative of a sym-polynomial
11 ;; Definite integral evaluation
12 #:evaluate-definite ; numeric result of ∫[a,b] p dx
13 #:evaluate-integral ; numeric result of a definite sym-integral
14 ;; Higher-order / iterated integrals
15 #:integrate-n ; n-th antiderivative
16 ;; sym-integral shell helpers
17 #:make-indefinite-integral ; wrap an integrand in an indefinite sym-integral
18 #:make-definite-integral ; wrap an integrand + bounds in a sym-integral
19 ;; Smoke test
20 #:run-smoke-test))
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)) ; => 1x^3 [3/(2+1) = 1]
11 (integrate-term (make-term -1 x 1)) ; => -1/2x^2
12 (integrate-term (make-term 5 x 0)) ; => 5x"
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 and the sum
4 rule across terms.
5
6 Examples:
7 poly = 3x² - x + 5
8 ∫poly = x³ - (1/2)x² + 5x
9
10 Returns a new sym-polynomial of degree (degree(POLY) + 1)."
11 (check-type poly sym-polynomial)
12 (let* ((var (polynomial-variable poly))
13 (new-terms (mapcar #'integrate-term (polynomial-terms poly))))
14 (if (null new-terms)
15 ;; integral of the zero polynomial is still zero
16 (zero-polynomial var :domain (polynomial-domain poly))
17 (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, and here too the domain is carried through to the result: integrating the complex zero polynomial returns a complex zero polynomial rather than silently resetting it to :real. 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) ; ∫₀¹ p dx
20 (evaluate-definite p 0 (make-constant 'pi :pi)) ; ∫₀^π p dx"
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. Given a bound it returns a number: plain numbers pass straight through, a sym-constant is resolved through constant-numeric-value, nil raises an error explaining that only definite integrals have a numeric value, and anything else raises an error naming the offending type. The two coerce calls in evaluate-definite then ensure the bounds are double-float before they reach 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.
Evaluating an Integral Object
The sym-integral shells built by make-definite-integral are convenient for display, but until now there was no way to get a number back out of one without unpacking it by hand. evaluate-integral closes that loop: it checks that the integral is definite, checks that the stored integrand is a polynomial, and hands the pieces to evaluate-definite.
1 (defun evaluate-integral (integral)
2 "Return the numeric value of the definite sym-integral INTEGRAL.
3
4 This is the bridge between the data layer's sym-integral shells and the
5 numeric evaluator: it unpacks the integrand and bounds and delegates to
6 EVALUATE-DEFINITE.
7
8 Example:
9 (evaluate-integral (make-definite-integral p 0 1)) ; => 5.5"
10 (check-type integral sym-integral)
11 (unless (integral-definite-p integral)
12 (error "INTEGRAL is indefinite; only definite integrals have a numeric value."))
13 (assert (polynomial-p (integral-integrand integral))
14 (integral)
15 "INTEGRAL's integrand must be a sym-polynomial, got ~s"
16 (type-of (integral-integrand integral)))
17 (evaluate-definite (integral-integrand integral)
18 (integral-lower integral)
19 (integral-upper integral)))
The two guards are what make the shells safe to pass around. An indefinite integral has no numeric value at all, so evaluate-integral says so explicitly instead of failing somewhere deep inside %resolve-bound; and because the integrand slot is intentionally untyped — recall that sym-integral accepts “any form” so that it can represent expressions the polynomial layer cannot yet handle — a non-polynomial integrand produces a clear assertion naming the actual type rather than a confusing check-type failure. In practice you would call evaluate-integral on a shell you had previously displayed with integral->string.
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 in indefinite form.
3 The stored integrand is POLY itself, not its antiderivative; call
4 INTEGRATE to take the antiderivative.
5
6 Rendering: ∫(expr) dx"
7 (check-type poly sym-polynomial)
8 (make-integral poly (polynomial-variable poly)))
9
10 (defun make-definite-integral (poly lower upper)
11 "Create a sym-integral wrapping POLY with explicit bounds LOWER and UPPER.
12 The stored integrand is POLY itself, not its antiderivative; call
13 EVALUATE-INTEGRAL to get the numeric value.
14 LOWER and UPPER may be real numbers or sym-constant objects.
15
16 Rendering: ∫[a,b](expr) dx"
17 (check-type poly sym-polynomial)
18 (make-integral poly (polynomial-variable poly) :lower lower :upper upper))
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
11 ;; ∫p dx = x³ - (1/2)x² + 5x
12 (ip (integrate p))
13
14 ;; ∫²p dx = (1/4)x⁴ - (1/6)x³ + (5/2)x²
15 (iip (integrate-n p 2))
16
17 ;; q = 6x + 2
18 (q (make-polynomial x (list (make-term 6 x 1)
19 (make-term 2 x 0))))
20
21 ;; ∫q dx = 3x² + 2x ← exact antiderivative
22 (iq (integrate q))
23
24 ;; Definite: ∫₀¹ p dx = F(1) - F(0)
25 ;; F(x) = x³ - (1/2)x² + 5x
26 ;; F(1) = 1 - 1/2 + 5 = 5.5
27 ;; F(0) = 0
28 ;; ⇒ 5.5
29 (def-01 (evaluate-definite p 0 1))
30
31 ;; Definite: ∫₀¹ q dx = [3x² + 2x]₀¹ = 3 + 2 = 5
32 (def-q-01 (evaluate-definite q 0 1))
33
34 ;; Definite: ∫₀^π p dx
35 (def-pi (evaluate-definite p 0 pi-c))
36
37 ;; sym-integral shells
38 (indef-shell (make-indefinite-integral p))
39 (def-shell (make-definite-integral p 0 1))
40 (pi-shell (make-definite-integral p 0 pi-c))
41
42 ;; evaluating a shell directly (no need to unpack it by hand)
43 (shell-01 (evaluate-integral def-shell)))
44
45 (format t "~%=== Integration Smoke Test ===~%~%")
46
47 (format t "p : ~a~%" (polynomial->string p))
48 (format t "∫p dx : ~a~%" (polynomial->string ip))
49 (format t "∫∫p dx dx : ~a~%" (polynomial->string iip))
50 (format t "~%")
51 (format t "q : ~a~%" (polynomial->string q))
52 (format t "∫q dx : ~a~%" (polynomial->string iq))
53 (format t "~%")
54 (format t "∫₀¹ p dx : ~a (expected 5.5)~%" def-01)
55 (format t "∫₀¹ q dx : ~a (expected 5.0)~%" def-q-01)
56 (format t "∫₀^π p dx : ~f~%" def-pi)
57 (format t "~%")
58 (format t "indef shell : ~a~%" (integral->string indef-shell))
59 (format t "def [0,1] : ~a~%" (integral->string def-shell))
60 (format t "def [0,pi] : ~a~%" (integral->string pi-shell))
61 (format t "~%")
62 (format t "evaluate shell : ~a (expected 5.5)~%" shell-01)
63 (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 (π), the string rendering of all three integral shell variants, and one evaluate-integral call that confirms a shell can be evaluated without being unpacked. The ~f directive for the π-bounded result prints a floating-point number whose exact digits vary a little between implementations; the value should be about 41.7794, which is π³ − π²/2 + 5π. 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 T
4 * (in-package #:symbolic-math/integ)
5 #<PACKAGE "SYMBOLIC-MATH/INTEG">
6 * (run-smoke-test)
7
8 === Integration Smoke Test ===
9
10 p : 3x^2 + -1x + 5
11 ∫p dx : 1x^3 + -1/2x^2 + 5x
12 ∫∫p dx dx : 1/4x^4 + -1/6x^3 + 5/2x^2
13
14 q : 6x + 2
15 ∫q dx : 3x^2 + 2x
16
17 ∫₀¹ p dx : 5.5d0 (expected 5.5)
18 ∫₀¹ q dx : 5.0d0 (expected 5.0)
19 ∫₀^π p dx : 41.77943774770411
20
21 indef shell : ∫(3x^2 + -1x + 5) dx
22 def [0,1] : ∫[0,1](3x^2 + -1x + 5) dx
23 def [0,pi] : ∫[0,π](3x^2 + -1x + 5) dx
24
25 evaluate shell : 5.5d0 (expected 5.5)
26
27 ==============================
28 NIL
29 *
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 — removes 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, 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.
Optional Practice Problems
Polynomial Multiplication Support: The data layer in data.lisp implements
polynomial-addandpolynomial-subtractbut does not support multiplication. Write a functionpolynomial-multiplythat multiplies two polynomialspandq. The algorithm should multiply each term ofpwith each term ofq(adding exponents and multiplying coefficients) and combine the resulting terms usingmake-polynomial.The Product Rule for Symbolic Differentiation: Currently, differentiation in differentiation.lisp only supports single polynomials by applying the power rule term-by-term. In real calculus, the derivative of a product of two functions is given by
(f·g)' = f'·g + f·g'. Write a functiondifferentiate-productthat takes twosym-polynomialstructures, computes their derivative using the product rule (leveragingdifferentiateand your newpolynomial-multiplyfunction), and returns the resulting simplified polynomial.Symbolic Integration by Substitution (Simple Linear Substitution): In integration.lisp, the reverse power rule is applied directly to monomials. Extend the integration engine to support simple u-substitution for linear compositions of the form
∫(ax + b)^n dx = (ax + b)^(n+1) / (a(n+1)). Represent the composition using a custom structure or by parsing list representation, and return the integrated polynomial.Finding Roots via Newton’s Method: Using the symbolic differentiation engine, we can compute exact derivatives at any point. Implement Newton’s method for finding roots of a polynomial:
x[n+1] = x[n] - f(x[n]) / f'(x[n]). Write a functionfind-polynomial-rootthat takes a polynomial, an initial guess, a tolerance, and a maximum iteration count. Usepolynomial-evaluateforf(x[n])andgradient-atfrom differentiation.lisp forf'(x[n]), and return the root.Definite Integral via Trapezoidal Rule (Numeric Verification): In integration.lisp, definite integrals are computed exactly using the Fundamental Theorem of Calculus. Write a numeric integration function,
evaluate-definite-trapezoidal, that estimates the definite integral of a polynomial over[a, b]using the trapezoidal rule withNsubintervals. Compare the numeric result against the exact result returned byevaluate-definiteto measure estimation error.Symbolic Math Parser for Lisp Expressions: The construction of polynomials like
3x² - x + 5is done using nested constructor calls. Write a macro or helper function,parse-math-expression, that takes a standard Lisp S-expression (e.g.(+ (* 3 (expt x 2)) (- x) 5)) and parses it into asym-polynomialobject. Ensure that variables, exponents, and constant terms are resolved correctly.