Symbolic Mathematics in Haskell

Symbolic computation manipulates mathematical expressions as data structures rather than as floating-point approximations. When you differentiate 3x^2 - x + 5 by hand, you apply the power rule to each term and write 6x - 1 then you are doing symbolic math. A numeric approach, by contrast, would estimate the derivative by evaluating the function at two nearby points and dividing. Symbolic computation gives you exact, algebraic results.

Haskell is a natural fit for symbolic math. Algebraic data types let you model mathematical expressions directly as trees of constructors. Pattern matching lets you write transformation rules that read almost exactly like the rules you learned in calculus. And because coefficients use Rational (exact fractions), differentiation and integration produce results with no floating-point drift.

Dear reader, you might wonder: why study symbolic math? In an age where everyone is talking about large language models there is room and a real need for deterministic technologies that stand on their own merits as well as augment sometime non-deterministic AI/LLM based systems. I hope you enjoy this material.

The library we will study in this chapter is located in the source-code/SymbolicMath directory, and provides symbolic differentiation and integration for single-variable polynomials. It was ported from the Common Lisp symbolic math library I wrote for my book Loving Common Lisp, or the Savvy Programmer’s Secret Weapon, and the translation from Lisp to Haskell is instructive: where the Lisp version uses lists and runtime type checking, the Haskell version uses statically-checked algebraic data types and precise function signatures.

The library is organized as a Cabal project with three exposed modules and a smoke-test executable:

1 SymbolicMath/
2 ├── symbolic-math.cabal
3 ├── src/SymbolicMath/
4 │   ├── Types.hs              -- Core data structures
5 │   ├── Differentiation.hs    -- Symbolic differentiation
6 │   └── Integration.hs        -- Symbolic integration
7 └── app/Main.hs               -- Smoke test runner

Build and run with:

1 cd source-code/SymbolicMath
2 cabal build
3 cabal run symbolic-math

The smoke test exercises every function in the library and prints expected versus actual results. If you are working through this chapter, keep a REPL open with cabal repl and experiment with the examples as you read.

The module SymbolicMath.Types defines five data types that together model everything we need for single-variable polynomial calculus. Let us walk through each one.

A Domain tags whether a variable or polynomial lives in the reals, the complex numbers, or the integers:

1 data Domain = Real | Complex | Integer_
2   deriving (Show, Eq)

The library uses Domain mostly for bookkeeping; the differentiation and integration rules do not inspect it. But encoding the domain in the type gives you the option to add domain-specific validation later. For example this code refuses to take a square root in Integer_ domain.

A Variable is a named symbol with a domain:

1 data Variable = Variable
2   { varName   :: String
3   , varDomain :: Domain
4   } deriving (Show)

Two variables are equal when both name and domain match:

1 makeVariable :: String -> Domain -> Variable
2 makeVariable name domain = Variable name domain

In the REPL:

1 ghci> let x = makeVariable "x" Real
2 ghci> let n = makeVariable "n" Integer_
3 ghci> variableEq x x
4 True
5 ghci> variableEq x (makeVariable "x" Complex)
6 False

A Constant represents a named mathematical constant like pi or e, or an exact rational number you want to treat symbolically:

 1 data ConstantValue
 2   = ExactValue Rational
 3   | SymbolicPi
 4   | SymbolicE
 5   deriving (Show, Eq)
 6 
 7 data Constant = Constant
 8   { constName  :: String
 9   , constValue :: ConstantValue
10   } deriving (Show, Eq)

The function constantNumericValue collapses a constant to a Double when you need a numeric answer:

1 ghci> let piConst = makeConstant "pi" SymbolicPi
2 ghci> let eConst  = makeConstant "e"  SymbolicE
3 ghci> let g       = makeConstant "g"  (ExactValue 9.80665)
4 ghci> constantNumericValue piConst
5 3.141592653589793
6 ghci> constantNumericValue g
7 9.80665

Constants are used most often as integral bounds such as writing int_0^pi reads more naturally than int_0^3.14159.

A Term is a single monomial: coefficient times variable raised to a non-negative integer exponent:

1 data Term = Term
2   { termCoefficient :: Rational
3   , termVariable    :: Variable
4   , termExponent    :: Int
5   } deriving (Show)

The constructor makeTerm enforces the non-negative exponent invariant:

1 makeTerm :: Rational -> Variable -> Int -> Term
2 makeTerm coeff var exp_
3   | exp_ < 0  = error "makeTerm: exponent must be non-negative"
4   | otherwise = Term coeff var exp_

A few examples:

1 ghci> makeTerm 3  x 2    -- 3x²
2 ghci> makeTerm (-1) x 1  -- -x
3 ghci> makeTerm 5  x 0    -- constant 5

The module provides termNegate (multiply coefficient by -1), termScale (multiply coefficient by a scalar), and termToString for display:

1 ghci> termToString (makeTerm 3 x 2)
2 "3x^2"
3 ghci> termToString (makeTerm (-1) x 1)
4 "-1x"
5 ghci> termToString (makeTerm 5 x 0)
6 "5"

A Polynomial is an ordered list of Term values in a single variable, with terms stored in descending exponent order:

1 data Polynomial = Polynomial
2   { polyVariable :: Variable
3   , polyTerms    :: [Term]
4   , polyDomain   :: Domain
5   } deriving (Show)

The smart constructor makePolynomial is where the real work happens. It takes a variable, a list of terms (possibly unsorted, possibly with duplicate exponents), and a domain. It combines like-exponent terms by summing their coefficients, drops any terms whose coefficient sums to zero, and sorts by descending exponent:

1 makePolynomial :: Variable -> [Term] -> Domain -> Polynomial
2 makePolynomial var terms domain =
3   let coeffMap = Map.fromListWith (+) [(termExponent t, termCoefficient t) | t <- terms]
4       pairs    = [(e, coeff) | (e, coeff) <- Map.toList coeffMap, coeff /= 0]
5       sorted   = sortOn (\(e, _) -> Down e) pairs
6       newTerms = [Term coeff var e | (e, coeff) <- sorted]
7   in Polynomial var newTerms domain

This uses Data.Map.Strict to group by exponent and sum coefficients which is a nice example of using the right data structure to simplify the logic. The Down newtype from Data.Ord gives us descending sort order.

Let us build a polynomial:

1 ghci> let x = makeVariable "x" Real
2 ghci> let p = makePolynomial x
3         [ makeTerm  3  x 2
4         , makeTerm (-1) x 1
5         , makeTerm  5  x 0 ] Real
6 ghci> polynomialToString p
7 "3x^2 + -1x + 5"

Since makePolynomial normalizes automatically, you can feed it messy input and still get a canonical result:

1 ghci> let messy = makePolynomial x
2         [ makeTerm 2 x 2, makeTerm 1 x 2, makeTerm 3 x 1
3         , makeTerm (-3) x 1, makeTerm 4 x 0 ] Real
4 ghci> polynomialToString messy
5 "3x^2 + 4"

The 2x^2 and 1x^2 combined to 3x^2; the 3x and -3x canceled; the constant 4 remained. This normalization guarantee means every polynomial has exactly one representation, which simplifies equality checking and debugging.

The module provides addition, subtraction, negation, and scalar multiplication. All return new polynomials and inputs are never mutated:

1 ghci> let q = makePolynomial x [makeTerm 1 x 1, makeTerm 2 x 0] Real
2 ghci> polynomialToString (polynomialAdd p q)
3 "3x^2 + 7"
4 ghci> polynomialToString (polynomialSubtract p q)
5 "3x^2 + -2x + 3"
6 ghci> polynomialToString (polynomialScale 2 p)
7 "6x^2 + -2x + 10"

Addition checks that both polynomials use the same variable and errors out if they do not so you cannot add a polynomial in x to a polynomial in y.

polynomialEvaluate substitutes a numeric value for the variable. The type signature uses Fractional a so it works with Rational, Double, or any fractional type:

1 polynomialEvaluate :: Fractional a => Polynomial -> a -> a
2 polynomialEvaluate (Polynomial _ terms _) x =
3   sum [fromRational c * x ^ n | Term c _ n <- terms]

1 ghci> polynomialEvaluate p (0 :: Rational)
2 5
3 ghci> polynomialEvaluate p (1 :: Rational)
4 7
5 ghci> polynomialEvaluate p (2 :: Rational)
6 15

polynomialDegree returns the highest exponent, or -1 for the zero polynomial (a convention borrowed from computer algebra systems). polynomialLeadingTerm returns the first term in the sorted list, or Nothing.

Three convenience constructors round out the API:

1 zeroPolynomial var        -- 0
2 constantPolynomial var 7  -- 7
3 identityPolynomial var    -- x

The final data type stores an unevaluated integral along with its bounds:

1 data SymIntegral = SymIntegral
2   { integrand   :: Polynomial
3   , intVariable :: Variable
4   , intLower    :: Maybe (Either Rational Constant)
5   , intUpper    :: Maybe (Either Rational Constant)
6   } deriving (Show, Eq)

Bounds use Maybe because an integral is either definite (both bounds present) or indefinite (both absent). The Either Rational Constant lets you write numeric bounds like 0 and 1 or symbolic bounds like pi and e. The constructor enforces the invariant that both bounds must be Just or both Nothing:

1 makeIntegral :: Polynomial -> Variable
2              -> Maybe (Either Rational Constant)
3              -> Maybe (Either Rational Constant)
4              -> SymIntegral

 1 ghci> -- indefinite: ∫(3x² - x + 5) dx
 2 ghci> let indef = makeIntegral p x Nothing Nothing
 3 ghci> -- definite: ∫₀¹ (3x² - x + 5) dx
 4 ghci> let def01 = makeIntegral p x (Just (Left 0)) (Just (Left 1))
 5 ghci> -- definite with symbolic bound: ∫₀^π
 6 ghci> let piInt = makeIntegral p x (Just (Left 0)) (Just (Right piConst))
 7 ghci> integralToString indef
 8 "∫(3x^2 + -1x + 5) dx"
 9 ghci> integralToString def01
10 "∫[0,1](3x^2 + -1x + 5) dx"

The module SymbolicMath.Differentiation implements the power rule and sum rule. The core function operates on a single term:

1 differentiateTerm :: Term -> Maybe Term
2 differentiateTerm (Term c v n)
3   | n == 0    = Nothing          -- derivative of a constant is 0
4   | otherwise = Just (Term (fromIntegral n * c) v (n - 1))

The Maybe return type handles the case where a term differentiates to zero so we want to drop those terms from the result polynomial rather than storing a zero-coefficient term. This is exactly the pattern the power rule describes:

d/dx (c · x^{n) = n · c · x}(n-1) for n ≥ 1

The polynomial-level differentiate applies differentiateTerm to each term via a local mapMaybe helper, then constructs a new normalized polynomial (or the zero polynomial if every term differentiated to zero):

1 differentiate :: Polynomial -> Polynomial
2 differentiate (Polynomial var terms domain) =
3   let diffTerms = mapMaybeT differentiateTerm terms
4   in if null diffTerms
5      then zeroPolynomial var
6      else makePolynomial var diffTerms domain

Higher-order derivatives use simple recursion:

1 differentiateN :: Polynomial -> Int -> Polynomial
2 differentiateN poly n
3   | n < 0     = error "differentiateN: n must be non-negative"
4   | n == 0    = poly
5   | otherwise = differentiateN (differentiate poly) (n - 1)

Two convenience functions bridge the symbolic and numeric worlds. gradientAt differentiates symbolically and then evaluates numerically at a point. criticalPointP checks whether the absolute value of the derivative at a point falls below a tolerance:

1 gradientAt :: Polynomial -> Double -> Double
2 gradientAt poly x = polynomialEvaluate (differentiate poly) x
3 
4 criticalPointP :: Polynomial -> Double -> Double -> Bool
5 criticalPointP poly x tolerance = abs (gradientAt poly x) < tolerance

Let us see the differentiation functions in action. We will use the polynomial p = 3x^2 - x + 5:

 1 ghci> let x = makeVariable "x" Real
 2 ghci> let p = makePolynomial x
 3         [ makeTerm 3  x 2
 4         , makeTerm (-1) x 1
 5         , makeTerm 5  x 0 ] Real
 6 ghci> let dp = differentiate p
 7 ghci> polynomialToString dp
 8 "6x + -1"
 9 ghci> let ddp = differentiate dp
10 ghci> polynomialToString ddp
11 "6"
12 ghci> let dddp = differentiate ddp
13 ghci> polynomialToString dddp
14 "0"
15 ghci> gradientAt p 0
16 -1.0
17 ghci> gradientAt p 1
18 5.0
19 ghci> criticalPointP p (1/6) 1e-9
20 True

The critical point check at x = 1/6 returns True because p’(x) = 6x - 1 equals zero at x = 1/6. This is the minimum of the parabola.

Let us also differentiate q = x^{4 - 2x}3 + x to see higher-degree behavior:

1 ghci> let q = makePolynomial x
2         [ makeTerm 1  x 4
3         , makeTerm (-2) x 3
4         , makeTerm 1  x 1 ] Real
5 ghci> polynomialToString (differentiate q)
6 "4x^3 + -6x^2 + 1"
7 ghci> polynomialToString (differentiateN q 2)
8 "12x^2 + -12x"

The module SymbolicMath.Integration implements the reverse power rule and the Fundamental Theorem of Calculus for definite integrals.

The reverse power rule on a single term is:

∫ c · x^{n dx = (c/(n+1)) · x}(n+1)

In Haskell:

1 integrateTerm :: Term -> Term
2 integrateTerm (Term c v n) =
3   let newExp  = n + 1
4       newCoef = c / fromIntegral newExp
5   in Term newCoef v newExp

Unlike differentiation, integration of a term never produces zero — even a constant term c integrates to c · x. So integrateTerm returns a Term directly rather than a Maybe Term.

The polynomial-level integrate maps integrateTerm over every term and normalizes:

1 integrate :: Polynomial -> Polynomial
2 integrate (Polynomial var terms domain) =
3   let newTerms = map integrateTerm terms
4   in if null newTerms
5      then zeroPolynomial var
6      else makePolynomial var newTerms domain

Note that integrate does not add a constant of integration +C. The antiderivative is returned as a polynomial, and you can add a constant term yourself if needed.

For definite integrals, evaluateDefinite uses the Fundamental Theorem of Calculus:

∫[a,b] f(x) dx = F(b) - F(a) where F = ∫ f

1 evaluateDefinite :: Polynomial
2                  -> Either Rational Constant
3                  -> Either Rational Constant
4                  -> Double
5 evaluateDefinite poly lower upper =
6   let f = integrate poly
7       a = resolveBoundToDouble lower
8       b = resolveBoundToDouble upper
9   in polynomialEvaluate f b - polynomialEvaluate f a

The helper resolveBoundToDouble converts either a rational or a symbolic constant to Double. This is the one place in the library where exact arithmetic gives way to floating-point and it happens at the very last step, after all symbolic manipulation is complete.

Higher-order integration uses the same recursive pattern as higher-order differentiation:

1 integrateN :: Polynomial -> Int -> Polynomial
2 integrateN poly n
3   | n < 0     = error "integrateN: n must be non-negative"
4   | n == 0    = poly
5   | otherwise = integrateN (integrate poly) (n - 1)

Two convenience functions wrap a polynomial’s antiderivative in a SymIntegral for display:

1 makeIndefiniteIntegral :: Polynomial -> SymIntegral
2 makeDefiniteIntegral :: Polynomial
3                      -> Either Rational Constant
4                      -> Either Rational Constant
5                      -> SymIntegral

Let us run through the integration examples using p = 3x^2 - x + 5:

 1 ghci> let x = makeVariable "x" Real
 2 ghci> let piC = makeConstant "pi" SymbolicPi
 3 ghci> let p = makePolynomial x
 4         [ makeTerm 3  x 2
 5         , makeTerm (-1) x 1
 6         , makeTerm 5  x 0 ] Real
 7 ghci> let ip = integrate p
 8 ghci> polynomialToString ip
 9 "1x^3 + -1/2x^2 + 5x"
10 ghci> evaluateDefinite p (Left 0) (Left 1)
11 5.5
12 ghci> evaluateDefinite p (Left 0) (Right piC)
13 41.779...

The indefinite integral of 3x^2 - x + 5 is x^{3 - (1/2)x}2 + 5x. The definite integral from 0 to 1 evaluates to 5.5 and you can verify this by hand: F(1) - F(0) = (1 - 0.5 + 5) - 0 = 5.5.

We can also integrate a simpler polynomial q = 6x + 2:

1 ghci> let q = makePolynomial x [makeTerm 6 x 1, makeTerm 2 x 0] Real
2 ghci> polynomialToString (integrate q)
3 "3x^2 + 2x"
4 ghci> evaluateDefinite q (Left 0) (Left 1)
5 5.0

The app/Main.hs file runs all three smoke tests in sequence:

 1 module Main where
 2 
 3 import qualified SymbolicMath.Types           as T
 4 import qualified SymbolicMath.Differentiation as D
 5 import qualified SymbolicMath.Integration     as I
 6 
 7 main :: IO ()
 8 main = do
 9   T.runSmokeTest
10   D.runSmokeTest
11   I.runSmokeTest

Running it produces:

 1 === Symbolic Math Data Layer Smoke Test ===
 2 
 3 p         : 3x^2 + -1x + 5
 4 q         : 1x + 2
 5 p + q     : 3x^2 + 7
 6 p - q     : 3x^2 + -2x + 3
 7 2 * p     : 6x^2 + -2x + 10
 8 degree(p) : 2
 9 p(0)      : 5
10 p(1)      : 7
11 p(2)      : 15
12 indef     : ∫(3x^2 + -1x + 5) dx
13 def       : ∫[0,1](3x^2 + -1x + 5) dx
14 pi-int    : ∫[0,pi](3x^2 + -1x + 5) dx
15 
16 === Differentiation Smoke Test ===
17 
18 p            : 3x^2 + -1x + 5
19 p'           : 6x + -1
20 p''          : 6
21 p'''         : 0
22 
23 q            : 1x^4 + -2x^3 + 1x
24 q'           : 4x^3 + -6x^2 + 1
25 q'' (via n=2): 12x^2 + -12x
26 
27 gradient-at(p, 0)      : -1.0  (expected -1)
28 gradient-at(p, 1)      : 5.0  (expected  5)
29 critical-point-p(p,1/6): True  (expected True)
30 critical-point-p(p,0)  : False (expected False)
31 
32 === Integration Smoke Test ===
33 
34 p              : 3x^2 + -1x + 5
35 ∫p dx          : 1x^3 + -1/2x^2 + 5x
36 ∫∫p dx dx      : 1/4x^4 + -1/6x^3 + 5/2x^2
37 
38 q              : 1x + 2
39 ∫q dx          : 3x^2 + 2x
40 
41 ∫₀¹  p dx      : 5.5  (expected 5.5)
42 ∫₀¹  q dx      : 5.0  (expected 5.0)
43 ∫₀^π p dx      : 41.779...

A few architectural choices are worth calling out because they make the library pleasant to work with and safe against common errors:

Exact arithmetic. Coefficients use Rational (a ratio of Integer values), so differentiation and integration produce exact fractions. Floating-point only enters when you call evaluateDefinite or polynomialEvaluate with a Double argument — and even then, the symbolic manipulation that precedes evaluation is exact.

Automatic normalization. makePolynomial combines like terms and sorts by descending exponent on every call. There is no way to construct a polynomial with duplicate-degree terms or zero-coefficient terms. This invariant means equality checking is structural and polynomialToString always produces readable output.

Immutability. Every arithmetic function returns a fresh polynomial. You never need to worry about a differentiation call silently modifying a polynomial you are still using elsewhere.

Bound flexibility. Integral bounds accept Rational numbers or Constant values, making it straightforward to express ∫[0,π] or ∫[0,e] without manually converting to floating-point.

Single-variable limitation. The library handles one variable at a time. Multi-variable polynomials and partial derivatives would require a different term representation (e.g., a map from tuples of variable-exponent pairs to coefficients). This is a natural extension if you want to deepen your understanding of the design.

We have built a small but complete symbolic math library in three modules. The key insight is that mathematical expressions are trees, and algebraic data types let you model those trees directly. Once you have the data types right, the transformation rules: the power rule, the sum rule, the reverse power rule, then you translate almost word-for-word into pattern-matching function definitions.

The library uses exact rational arithmetic throughout, normalizes polynomials automatically, and cleanly separates the symbolic manipulation layer from the numeric evaluation layer. The same design pattern of modeling your domain as algebraic data types, then writing pure functions that transform those types. This pattern applies to compilers, theorem provers, equation solvers, and many other domains beyond calculus.

I encourage you to experiment with the code in a REPL as you work through the practice problems below. The best way to internalize these ideas is to extend the library yourself.

1. Extend the Polynomial data type and arithmetic functions to support polynomials in two variables (e.g., 2x^2 y + 3xy - y). You will need to change the term representation so that a Term can hold an exponent for each variable. Then implement partialDerivativeX and partialDerivativeY functions that differentiate with respect to one variable while treating the other as a constant.

2. Add a simplify function to SymbolicMath.Types that improves the display of polynomials. Currently polynomialToString renders -1x + 5 with a leading coefficient of -1 and a trailing + before negative terms. Write simplify to produce cleaner output: render 1x as just x, omit the coefficient when it is 1 (except for constant terms), and replace + - with - between terms. For example, 3x^2 + -1x + 5 should become 3x^2 - x + 5.