The Functional Programmer
0.00
recursion26tagscollectionsarchivewikifeed

Building a CAS in OCaml Part 2

2026-04-13 · 10 min

ocamlcasalgorithms

Part 1 built integers that do not overflow. This builds exact rationals on them, then the expression tree, then the piece that makes the whole thing a CAS rather than a data structure: the automatic simplifier, whose job is to give every mathematically equal expression the same shape.

§ 01

Rationals

One invariant, maintained by a single normalizing constructor: positive denominator, sign on the numerator, coprime pair. Zero is exactly 0/1. That is what lets structural equality of the two fields mean mathematical equality, which the expression layer depends on when it collects like terms.

(* lib/rational.ml *)
 
type t = {
  num : Bigint.t;
  den : Bigint.t; (* always > 0, coprime with num *)
}
 
let normalize (num : Bigint.t) (den : Bigint.t) : t =
  if Bigint.is_zero den then raise Division_by_zero;
  if Bigint.is_zero num then { num = Bigint.zero; den = Bigint.one }
  else begin
    let s = Bigint.sign num * Bigint.sign den in
    let n = Bigint.abs num and d = Bigint.abs den in
    let g = Bigint.gcd n d in
    let n = Bigint.div n g and d = Bigint.div d g in
    { num = (if s < 0 then Bigint.neg n else n); den = d }
  end

Arithmetic is the schoolbook formulas, each routed back through normalize so no unreduced fraction ever escapes.

let add (a : t) (b : t) : t =
  normalize
    (Bigint.add (Bigint.mul a.num b.den) (Bigint.mul b.num a.den))
    (Bigint.mul a.den b.den)
 
let sub (a : t) (b : t) : t = add a (neg b)
let mul (a : t) (b : t) : t =
  normalize (Bigint.mul a.num b.num) (Bigint.mul a.den b.den)
 
let div (a : t) (b : t) : t =
  if is_zero b then raise Division_by_zero
  else normalize (Bigint.mul a.num b.den) (Bigint.mul a.den b.num)

Comparison cross-multiplies, which is valid precisely because both denominators are positive - the invariant paying off again.

let compare (a : t) (b : t) : int =
  Bigint.compare (Bigint.mul a.num b.den) (Bigint.mul b.num a.den)
 
let equal (a : t) (b : t) : bool =
  Bigint.equal a.num b.num && Bigint.equal a.den b.den
 
let pow (r : t) (n : int) : t =
  if n = min_int then invalid_arg "Rational.pow: exponent too large";
  if n >= 0 then normalize (Bigint.pow r.num n) (Bigint.pow r.den n)
  else begin
    if is_zero r then raise Division_by_zero;
    let m = -n in
    normalize (Bigint.pow r.den m) (Bigint.pow r.num m)
  end

The test that says why any of this is worth it: a third plus a third plus a third is exactly one, and the tenths add exactly, neither of which holds in binary floating point.

let test_exactness () =
  let third = frac 1 3 in
  let sum = Rational.add (Rational.add third third) third in
  check_str "1/3 three times is exactly 1" "1" (s sum);
  check "tenths add exactly"
    (Rational.equal (Rational.add (frac 1 10) (frac 2 10)) (frac 3 10));
  let h = ref Rational.zero in
  for i = 1 to 10 do
    h := Rational.add !h (frac 1 i)
  done;
  check_str "H_10" "7381/2520" (s !h)
§ 02

The expression tree

Sums and products are n-ary, not binary. That is the single most consequential design choice here: a + b + c is one node with three operands rather than two nested nodes, so flattening and sorting are list operations and there is no left-versus-right associativity to canonicalize away.

(* lib/expr.ml *)
 
exception Undefined of string
 
type t =
  | Num of Rational.t
  | Sym of string
  | Add of t list          (* n-ary, flattened, sorted, >= 2 operands *)
  | Mul of t list          (* n-ary, flattened, sorted, >= 2 operands *)
  | Pow of t * t
  | Fun of string * t list (* argument order is significant: never sorted *)

A canonical order. Any total, deterministic order produces a canonical form - this one is picked to be cheap, and to rank numbers first so a product's coefficient is always its head, which the term-splitting below relies on.

let rank (e : t) : int =
  match e with
  | Num _ -> 0
  | Sym _ -> 1
  | Fun _ -> 2
  | Pow _ -> 3
  | Mul _ -> 4
  | Add _ -> 5
 
let rec compare_expr (a : t) (b : t) : int =
  match (a, b) with
  | Num x, Num y -> Rational.compare x y
  | Sym x, Sym y -> String.compare x y
  | Pow (b1, e1), Pow (b2, e2) ->
    let c = compare_expr b1 b2 in
    if c <> 0 then c else compare_expr e1 e2
  | Fun (f1, a1), Fun (f2, a2) ->
    let c = String.compare f1 f2 in
    if c <> 0 then c else compare_list a1 a2
  | Add x, Add y -> compare_list x y
  | Mul x, Mul y -> compare_list x y
  | _ -> compare (rank a) (rank b)
 
and compare_list (xs : t list) (ys : t list) : int =
  match (xs, ys) with
  | [], [] -> 0
  | [], _ -> -1
  | _, [] -> 1
  | x :: xs, y :: ys ->
    let c = compare_expr x y in
    if c <> 0 then c else compare_list xs ys
§ 03

Splitting terms and factors

Collecting like terms needs a way to ask what two terms have in common. A term of a sum splits into a rational coefficient and the rest; a factor of a product splits into a base and an exponent. Both splits are cheap because operands are already sorted.

3*x*y and -2*x*y split to the same term part, so they can be added. Numbers rank first, so any coefficient is the head of the product - no scan required.

let split_coeff (e : t) : Rational.t * t =
  match e with
  | Mul (Num c :: rest) ->
    let term = match rest with [] -> num_one | [ x ] -> x | xs -> Mul xs in
    (c, term)
  | _ -> (Rational.one, e)
 
let split_power (e : t) : t * t =
  match e with Pow (b, x) -> (b, x) | _ -> (e, num_one)
 
let rebuild_term (c : Rational.t) (term : t) : t =
  if Rational.is_one c then term
  else if is_num_one term then Num c
  else match term with Mul xs -> Mul (Num c :: xs) | _ -> Mul [ Num c; term ]
§ 04

Simplifying a sum

Flatten nested sums, fold every numeric operand into one constant, split the rest into coefficient-and-term pairs, sort so equal term parts become adjacent, merge them in one pass, drop anything whose coefficient cancelled to zero, and unwrap a one-element sum.

and simplify_sum (terms : t list) : t =
  let terms = List.concat_map (function Add xs -> xs | x -> [ x ]) terms in
  let const, pairs =
    List.fold_left
      (fun (c, acc) x ->
        match x with
        | Num r -> (Rational.add c r, acc)
        | _ -> (c, split_coeff x :: acc))
      (Rational.zero, []) terms
  in
  let sorted =
    List.sort (fun (_, t1) (_, t2) -> compare_expr t1 t2) (List.rev pairs)
  in
  let collected = collect_sum sorted in
  let others =
    List.filter_map
      (fun (c, term) ->
        if Rational.is_zero c then None else Some (rebuild_term c term))
      collected
  in
  let all = if Rational.is_zero const then others else Num const :: others in
  match all with
  | [] -> num_zero
  | [ x ] -> x
  | xs -> Add (List.sort compare_expr xs)

Sorting is what makes collection a single linear pass rather than a quadratic search: equal term parts are necessarily adjacent.

and collect_sum (pairs : (Rational.t * t) list) : (Rational.t * t) list =
  match pairs with
  | (c1, t1) :: (c2, t2) :: rest when compare_expr t1 t2 = 0 ->
    collect_sum ((Rational.add c1 c2, t1) :: rest)
  | x :: rest -> x :: collect_sum rest
  | [] -> []
§ 05

Simplifying a product

The same shape, with two differences: a single zero factor short-circuits the whole product, and recombining exponents can collapse a factor back to a plain number, which then has to rejoin the constant rather than sit inside the product as a stray Num.

and simplify_product (factors : t list) : t =
  let factors = List.concat_map (function Mul xs -> xs | x -> [ x ]) factors in
  if List.exists is_num_zero factors then num_zero
  else begin
    let const, pairs =
      List.fold_left
        (fun (c, acc) x ->
          match x with
          | Num r -> (Rational.mul c r, acc)
          | _ -> (c, split_power x :: acc))
        (Rational.one, []) factors
    in
    let sorted =
      List.sort (fun (b1, _) (b2, _) -> compare_expr b1 b2) (List.rev pairs)
    in
    let collected = collect_product sorted in
    (* Recombining exponents can collapse a factor to a number -
       2^x * 2^(1-x) becomes 2^1 - so fold those back into the
       constant rather than leaving a Num inside the product. *)
    let const, others =
      List.fold_left
        (fun (c, acc) (b, x) ->
          match simplify_power b x with
          | Num r -> (Rational.mul c r, acc)
          | p -> (c, p :: acc))
        (const, []) collected
    in
    let others = List.rev others in
    if Rational.is_zero const then num_zero
    else begin
      let all = if Rational.is_one const then others else Num const :: others in
      match all with
      | [] -> num_one
      | [ x ] -> x
      | xs -> Mul (List.sort compare_expr xs)
    end
  end
 
and collect_product (pairs : (t * t) list) : (t * t) list =
  match pairs with
  | (b1, e1) :: (b2, e2) :: rest when compare_expr b1 b2 = 0 ->
    collect_product ((b1, simplify_sum [ e1; e2 ]) :: rest)
  | x :: rest -> x :: collect_product rest
  | [] -> []

That fold-back path is not hypothetical, and the test suite pins it.

2^x * 2^y        ->  2^(x + y)      exponents combine, stays symbolic
2^x * 2^(1-x)    ->  2              exponent collapses to 1, so the
                                    factor becomes a number and has to
                                    merge into the constant
§ 06

Simplifying a power, and the rules that are wrong

Two of the obvious power identities are false in general, and a CAS that applies them anyway will confidently produce wrong answers. Both are only valid when the outer exponent is an integer.

The counterexample worth keeping in mind: distributing a one-half exponent over a square turns 1 into -1.

((-1)^2)^(1/2)  =  1^(1/2)   =  1
(-1)^(2 * 1/2)  =  (-1)^1    =  -1
 
so (u^v)^w = u^(v*w) needs w to be an integer, and the same
restriction applies to (u*v)^w = u^w * v^w

The whole rule set, with the undefined cases raised rather than papered over.

and simplify_power (b : t) (x : t) : t =
  match (b, x) with
  (* A rational to an integer power is computed outright. *)
  | Num bb, Num xx when Rational.is_integer xx -> (
    match Rational.to_int_opt xx with
    | None -> Pow (b, x) (* exponent too large to materialize *)
    | Some n ->
      if Rational.is_zero bb && n = 0 then raise (Undefined "0^0")
      else if Rational.is_zero bb && n < 0 then raise (Undefined "0^negative")
      else Num (Rational.pow bb n))
  | _ ->
    if is_num_zero x then
      if is_num_zero b then raise (Undefined "0^0") else num_one
    else if is_num_one x then b
    else if is_num_zero b then
      (* 0^y for symbolic y: cannot decide the sign of y, so leave it. *)
      Pow (b, x)
    else if is_num_one b then num_one
    else (
      match b with
      | Pow (u, v) when is_integer_expr x ->
        simplify_power u (simplify_product [ v; x ])
      | Mul factors when is_integer_expr x ->
        simplify_product (List.map (fun f -> simplify_power f x) factors)
      | _ -> Pow (b, x))

The top-level recursion that ties the three together - simplify children first, then apply the rule for this node.

let rec simplify (e : t) : t =
  match e with
  | Num _ | Sym _ -> e
  | Add xs -> simplify_sum (List.map simplify xs)
  | Mul xs -> simplify_product (List.map simplify xs)
  | Pow (b, x) -> simplify_power (simplify b) (simplify x)
  | Fun (f, args) -> Fun (f, List.map simplify args)
§ 07

What canonical buys

Equality becomes one structural comparison of two simplified trees - no search, no rewriting to a common form at comparison time.

let equal (a : t) (b : t) : bool = compare_expr (simplify a) (simplify b) = 0
 
let add a b = simplify (Add [ a; b ])
let sub a b = simplify (Add [ a; Mul [ Num Rational.minus_one; b ] ])
let mul a b = simplify (Mul [ a; b ])
let div a b = simplify (Mul [ a; Pow (b, Num Rational.minus_one) ])
let neg a = simplify (Mul [ Num Rational.minus_one; a ])

Subtraction and division are not constructors. They are a product with -1 and a power of -1, which is why x - x and x / x collapse for free rather than needing their own cancellation rules.

x - x   is   Add [x; Mul [-1; x]]        -> coefficients 1 and -1 sum to 0  -> 0
x / x   is   Mul [x; Pow (x, -1)]        -> exponents 1 and -1 sum to 0     -> 1
x^2 / x is   Mul [Pow (x,2); Pow (x,-1)] -> exponents 2 and -1 sum to 1     -> x

The real output of the driver, which is the whole simplifier in one screen.

$ ./_out/main simplify

Verbatim.

x + x            ->  2*x
2*x + 3*x        ->  5*x
x * x * x        ->  x^3
x^2 * x^3        ->  x^5
x - x            ->  0
x / x            ->  1
(x*y)^2          ->  x^2*y^2
(x^2)^3          ->  x^6
1/2 + 1/3        ->  5/6
2^x * 2^(1-x)    ->  2
x + y + x        ->  y + 2*x
0 * x + 1 * y    ->  y

Two expressions built differently, landing on the same tree - which is what the canonical form is for.

$ ./_out/main demo

Verbatim. Note the printed order: the canonical order ranks a bare symbol before a product, so y comes before 5*x. It is arbitrary but total, which is all canonicalization requires - a nicer display order is a separate concern from deciding equality.

a           = y + 5*x
b           = y + 5*x
a equals b  = true
a tree      = (+ y (* 5 x))
a at x=2,y=1 = 11
§ 08

Printing

Precedence-driven, with two touches that matter for readability: a negative leading coefficient prints as a minus sign rather than a parenthesized factor, and a sum whose term already starts with a minus joins with ' - ' instead of ' + -'.

let prec (e : t) : int =
  match e with
  | Num r -> if Rational.sign r < 0 || not (Rational.is_integer r) then 2 else 4
  | Sym _ | Fun _ -> 4
  | Pow _ -> 3
  | Mul _ -> 2
  | Add _ -> 1
 
and raw (e : t) : string =
  match e with
  | Mul (Num c :: rest) when Rational.sign c < 0 ->
    let body = String.concat "*" (List.map (to_str 2) rest) in
    if Rational.equal c Rational.minus_one then "-" ^ body
    else "-" ^ Rational.to_string (Rational.abs c) ^ "*" ^ body
  | Add (first :: rest) ->
    let buf = Buffer.create 32 in
    Buffer.add_string buf (to_str 1 first);
    List.iter
      (fun term ->
        let s = to_str 1 term in
        if String.length s > 0 && s.[0] = '-' then begin
          Buffer.add_string buf " - ";
          Buffer.add_string buf (String.sub s 1 (String.length s - 1))
        end
        else begin
          Buffer.add_string buf " + ";
          Buffer.add_string buf s
        end)
      rest;
    Buffer.contents buf
  | _ -> (* ... *) ""
§ 09

Testing a canonical form

Example-based tests catch the rules you thought of. Three properties catch the ones you did not, and they are what actually found the bugs while writing this.

Idempotence. If any rule leaves behind output that the same rules would rewrite again, the second pass differs from the first - which is exactly the definition of a non-canonical result.

let test_idempotent () =
  Random.init 20260413;
  for _ = 1 to 2000 do
    let e = random_expr 3 in
    match simplify e with
    | exception Undefined _ -> ()   (* 0^0 is legitimately undefined *)
    | exception Division_by_zero -> ()
    | s ->
      let s2 = simplify s in
      if compare_expr s s2 <> 0 then begin
        incr failures;
        Printf.printf "FAIL: not idempotent\n  once  %s\n  twice %s\n"
          (to_sexp s) (to_sexp s2)
      end
  done

Order invariance. Reordering the operands of a sum or product must not change the simplified result - if it does, the output is not canonical, whatever else it is.

    match (simplify (Add [ a; b; c ]), simplify (Add [ c; a; b ])) with
    | exception Undefined _ -> ()
    | exception Division_by_zero -> ()
    | s1, s2 ->
      if compare_expr s1 s2 <> 0 then begin
        incr failures;
        Printf.printf "FAIL: sum not order-invariant\n  %s\n  %s\n"
          (to_sexp s1) (to_sexp s2)
      end

Value preservation, the one that catches a rule that is canonical but wrong: substitute numbers for every symbol in both the original and the simplified form, and the two must evaluate to the same rational.

let test_simplify_preserves_value () =
  Random.init 5150;
  for _ = 1 to 500 do
    let e = random_expr 3 in
    let xv = int (1 + Random.int 4) and yv = int (1 + Random.int 4) in
    let evaluate g = to_rational (substitute "y" yv (substitute "x" xv g)) in
    match (evaluate e, evaluate (simplify e)) with
    | exception Undefined _ -> ()
    | exception Division_by_zero -> ()
    | Some a, Some b ->
      if not (Rational.equal a b) then begin
        incr failures;
        Printf.printf "FAIL: simplification changed the value\n"
      end
    | _ -> ()
  done
§ 10

What is deliberately not here

Automatic simplification is the set of rewrites applied to everything, always. Expansion and factoring are not in it, and not by oversight: either can make an expression dramatically larger, so both are operations a user asks for by name. (x + y)^2 stays as written.

in scope, applied always:    constant folding, flattening, like terms,
                             like factors, identity elimination, sorting
 
out of scope, on request:    expand, factor, collect, together, cancel,
                             trigonometric and logarithmic identities
 
(x + y)^2  stays  (x + y)^2  until Part 3 gives you expand
§ 11

Download

The snapshot at the end of this part: bigint.ml, rational.ml, expr.ml, test_rational.ml, test_expr.ml, test_all.ml, main.ml. bigint.ml is unchanged from Part 1, so a diff of the two directories is exactly what this post added.

referenced by

related

← older

Building a CAS in OCaml Part 1

2026-04-06

newer →

Building a CAS in OCaml Part 3

2026-04-20