Building a CAS in OCaml Part 5

2026-05-04 · 14 min

Part 4 can factor and cancel, and it still cannot differentiate sin(x). This part adds the calculus: a derivative table, Laurent series, limits read off a valuation, and symbolic integration, which is the first operation in the series that is not a decision procedure and has to admit it.

§ 01

What Part 3 refused

Two refusals, both lifted here. The first by writing u^v as exp(v log u), the second by a table.

Part 3            diff sin(x), x   ->  no derivative rule for sin
                  diff x^y, x      ->  power with a symbolic exponent
 
Part 5            diff sin(x^2), x ->  2*x*cos(x^2)
                  diff x^x, x      ->  x^x*(1 + log(x))

The table. Every rule is written with constructors Expr already has, so nothing introduces a function the simplifier has never seen.

let table (f : string) (args : Expr.t list) : Expr.t =
  let u = match args with [ a ] -> a | _ -> raise (Not_differentiable (f ^ ": arity")) in
  let inv e = Expr.Pow (e, Expr.int (-1)) in
  let sq e = Expr.Pow (e, Expr.int 2) in
  match f with
  | "exp" -> Expr.Fun ("exp", [ u ])
  | "log" -> inv u
  | "sin" -> Expr.Fun ("cos", [ u ])
  | "cos" -> Expr.Mul [ Expr.int (-1); Expr.Fun ("sin", [ u ]) ]
  | "tan" -> Expr.Add [ Expr.int 1; sq (Expr.Fun ("tan", [ u ])) ]
  | "sinh" -> Expr.Fun ("cosh", [ u ])
  | "cosh" -> Expr.Fun ("sinh", [ u ])
  | "tanh" -> Expr.Add [ Expr.int 1; Expr.Mul [ Expr.int (-1); sq (Expr.Fun ("tanh", [ u ])) ] ]
  (* 1 / sqrt(1 - u^2) *)
  | "asin" -> Expr.Pow (Expr.Add [ Expr.int 1; Expr.Mul [ Expr.int (-1); sq u ] ], Expr.rat (-1) 2)
  | "acos" ->
    Expr.Mul
      [ Expr.int (-1);
        Expr.Pow (Expr.Add [ Expr.int 1; Expr.Mul [ Expr.int (-1); sq u ] ], Expr.rat (-1) 2) ]
  | "atan" -> inv (Expr.Add [ Expr.int 1; sq u ])
  | "sqrt" -> Expr.Mul [ Expr.rat 1 2; Expr.Pow (u, Expr.rat (-1) 2) ]
  | _ -> raise (Not_differentiable ("no derivative rule for " ^ f))

§ 02

The general power rule

A constant exponent is the power rule. A symbolic one needs u^v = exp(v log u), and differentiating that gives

\frac{d}{d x} u^{v} = u^{v} (v^{'} lo g u + \frac{v u ^{'}}{u})

Reaching for it only when the exponent really is symbolic keeps the common case free of logarithms.

  (* A constant exponent is the power rule and nothing more. *)
  | Expr.Pow (b, (Expr.Num n as ne)) ->
    Expr.Mul [ ne; Expr.Pow (b, Expr.Num (Rational.sub n Rational.one)); raw b v ]
  (* u^v = exp(v log u), so the general rule is
     u^v * (v' log u + v u' / u). Reaching for it only when the exponent
     really is symbolic keeps the common case free of logarithms. *)
  | Expr.Pow (b, x) ->
    let du = raw b v and dx = raw x v in
    Expr.Mul
      [ Expr.Pow (b, x);
        Expr.Add
          [ Expr.Mul [ dx; Expr.Fun ("log", [ b ]) ];
            Expr.Mul [ x; du; Expr.Pow (b, Expr.int (-1)) ] ] ]

Verbatim.

> diff x^x, x
x^x*(1 + log(x))
> diff a^x, x
a^x*log(a)
> diff exp(x)*log(x), x
exp(x)*log(x) + exp(x)*x^(-1)
> diff x^5, x, 3
60*x^2

§ 03

Laurent series

A series here is a valuation and a coefficient array: the value is sum c_i x^(v+i). The valuation is what makes it Laurent rather than Taylor, and it is the whole reason limits work. sin(x)/x is not a power series at all, but it is a Laurent series with valuation 0, and its constant term is the limit.

Coefficients are rationals, so a coefficient that is zero is zero rather than 1e-17.

type t = {
  v : int;                 (* valuation: the exponent of the first term *)
  c : Rational.t array;    (* c.(0) is the coefficient of x^v *)
}
 
let order (s : t) : int = Array.length s.c
 
(* Drop leading zero coefficients, raising the valuation to match. A
   series that is entirely zero to the known order has no meaningful
   valuation, so it is reported as zero at the truncation point. *)
let normalize (s : t) : t =
  let n = Array.length s.c in
  let i = ref 0 in
  while !i < n && Rational.is_zero s.c.(!i) do incr i done;
  if !i = 0 then s
  else if !i = n then { v = s.v + n; c = [||] }
  else { v = s.v + !i; c = Array.sub s.c !i (n - !i) }

Two series are aligned over a common valuation before anything binary happens, and the known length shrinks to whatever both sides guarantee.

(* Re-express two series over a common valuation and length. *)
let align (a : t) (b : t) : int * int * Rational.t array * Rational.t array =
  let v = min a.v b.v in
  let n = min (a.v + order a) (b.v + order b) - v in
  let n = max n 0 in
  (v, n,
   Array.init n (fun i -> coeff a (v + i)),
   Array.init n (fun i -> coeff b (v + i)))
 
let add (a : t) (b : t) : t =
  let v, n, x, y = align a b in
  normalize { v; c = Array.init n (fun i -> Rational.add x.(i) y.(i)) }

The reciprocal falls out of a * a^-1 = 1 as a recurrence. The leading coefficient has to be nonzero, which normalize guarantees for anything that is not zero to the known order.

(* Reciprocal, by the recurrence that falls out of a * a^-1 = 1. The
   leading coefficient must be nonzero, which normalize guarantees for
   anything that is not zero to the known order. *)
let inv (a : t) : t =
  if is_zero a then raise Division_by_zero;
  let a = normalize a in
  let n = order a in
  let c = Array.make n Rational.zero in
  c.(0) <- Rational.inv a.c.(0);
  for i = 1 to n - 1 do
    let acc = ref Rational.zero in
    for j = 1 to i do
      acc := Rational.add !acc (Rational.mul a.c.(j) c.(i - j))
    done;
    c.(i) <- Rational.neg (Rational.mul c.(0) !acc)
  done;
  { v = -a.v; c }

§ 04

Elementary functions, by recurrence

Defining exp by its Taylor coefficients would be a table; defining it by its differential equation makes composition free. For a series u with positive valuation, y = exp(u) satisfies y' = u' y, and matching coefficients gives

n y_{n} = k = 1 \sum n k u_{k} y_{n - k}

The positive-valuation requirement is not a limitation to work around: exp of a nonzero constant is not rational, so there is nothing to return.

let require_positive name (a : t) =
  let a = normalize a in
  if (not (is_zero a)) && a.v < 1 then
    raise (Cannot_expand (name ^ ": needs a series with zero constant term"));
  a
 
let exp_s (a : t) : t =
  let a = require_positive "exp" a in
  let n = order a in
  if n = 0 then const Rational.one 0
  else begin
    let y = Array.make n Rational.zero in
    y.(0) <- Rational.one;
    let da = diff a in
    (* n y_n = sum_{k=1..n} k a_k y_{n-k} *)
    for i = 1 to n - 1 do
      let acc = ref Rational.zero in
      for k = 1 to i do
        acc := Rational.add !acc (Rational.mul (Rational.mul (Rational.of_int k) (coeff a k)) y.(i - k))
      done;
      ignore da;
      y.(i) <- Rational.div !acc (Rational.of_int i)
    done;
    { v = 0; c = y }
  end

Checked against the textbook expansions. Compiled with ocamlopt 4.14.1.

exp(x)     = 1 + x + 1/2*x^2 + 1/6*x^3 + 1/24*x^4 + 1/120*x^5 + 1/720*x^6 + 1/5040*x^7 + O(x^8)
sin(x)     = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + O(x^9)
cos(x)     = 1 - 1/2*x^2 + 1/24*x^4 - 1/720*x^6 + O(x^8)
log(1+x)   = x - 1/2*x^2 + 1/3*x^3 - 1/4*x^4 + 1/5*x^5 - 1/6*x^6 + 1/7*x^7 - 1/8*x^8 + O(x^9)
sqrt(1+x)  = 1 + 1/2*x - 1/8*x^2 + 1/16*x^3 - 5/128*x^4 + 7/256*x^5 - 21/1024*x^6 + 33/2048*x^7 + O(x^8)
tan-ish    = x + 1/3*x^3 + 2/15*x^5 + 17/315*x^7 + O(x^9)
exp(sin x) = 1 + x + 1/2*x^2 - 1/8*x^4 - 1/15*x^5 - 1/240*x^6 + 1/90*x^7 + O(x^8)

The identities are the better test, because a reference value and an implementation can be wrong together in a way an identity cannot.

(* The identities the expansions have to satisfy, which catch an error
   the reference values would not if both were wrong the same way. *)
let test_identities () =
  let sin_ = Series.sin_s x and cos_ = Series.cos_s x and exp_ = Series.exp_s x in
  check "sin^2 + cos^2 = 1"
    (Series.is_zero (Series.sub (Series.add (Series.mul sin_ sin_) (Series.mul cos_ cos_)) one));
  check "d/dx sin = cos" (Series.is_zero (Series.sub (Series.diff sin_) cos_));
  check "d/dx cos = -sin" (Series.is_zero (Series.add (Series.diff cos_) sin_));
  check "d/dx exp = exp" (Series.is_zero (Series.sub (Series.diff exp_) exp_));
  check "exp(x)*exp(-x) = 1"
    (Series.is_zero (Series.sub (Series.mul exp_ (Series.exp_s (Series.neg x))) one));
  check "log(1+x) then exp gives 1+x"
    (Series.is_zero (Series.sub (Series.exp_s (Series.log1p_s x)) (Series.add one x)));
  check "sqrt(1+x) squared is 1+x"
    (let q = Series.binomial_s x (r 1 2) n in
     Series.is_zero (Series.sub (Series.mul q q) (Series.add one x)));
  check "integrating the derivative restores the series"
    (Series.is_zero (Series.sub (Series.integrate (Series.diff sin_)) sin_))

§ 05

Where the known order goes

Accumulating into a zero series at valuation 0 silently throws away terms: alignment takes the minimum valuation, so adding a series at valuation 1 to zero at valuation 0 shortens the result by one. The first version of sin did exactly that, and sin(x)/x came back one term weaker than it should have been.

Starting from the first term rather than from zero is the whole fix.

let sin_s (a : t) : t =
  let a = require_positive "sin" a in
  let n = order a in
  if is_zero a then a
  else begin
    (* Accumulate starting from the first term rather than from a zero
       series at valuation 0: adding to zero would align the result down
       to valuation 0 and throw away known terms at the top. *)
    let a2 = mul a a in
    let acc = ref a and term = ref a and k = ref 1 in
    let fact = ref (Rational.of_int 6) in
    while 2 * !k + 1 <= n && not (is_zero !term) do
      term := mul !term a2;
      let s = scale !term (Rational.inv !fact) in
      acc := (if !k mod 2 = 1 then sub !acc s else add !acc s);
      incr k;
      fact := Rational.mul !fact (Rational.of_int ((2 * !k) * ((2 * !k) + 1)))
    done;
    !acc
  end

Before and after, on the same input.

before   sin(x)   = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + O(x^8)
         sin(x)/x = 1 - 1/6*x^2 + 1/120*x^4 + O(x^7)
 
after    sin(x)   = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + O(x^9)
         sin(x)/x = 1 - 1/6*x^2 + 1/120*x^4 - 1/5040*x^6 + O(x^8)

§ 06

Expanding an expression

The bridge from Expr to Series. A fractional power only has a Laurent expansion when the base is 1 plus something vanishing, which is why x^(1/2) is refused rather than approximated.

let rec expand (e : Expr.t) (v : string) (n : int) : Series.t =
  match e with
  | Expr.Num r -> Series.const r n
  | Expr.Sym s -> if s = v then Series.ident n else raise (Cannot_expand ("free symbol " ^ s))
  | Expr.Add xs -> List.fold_left (fun acc x -> Series.add acc (expand x v n)) (Series.zero n) xs
  | Expr.Mul xs ->
    List.fold_left (fun acc x -> Series.mul acc (expand x v n)) (Series.const Rational.one n) xs
  | Expr.Pow (b, Expr.Num r) when Rational.is_integer r -> (
    match Rational.to_int_opt r with
    | Some k -> Series.pow_int (expand b v n) k
    | None -> raise (Cannot_expand "exponent too large"))
  (* A fractional power only has a Laurent expansion when the base is
     1 + (something vanishing); x^(1/2) is not a Laurent series at all. *)
  | Expr.Pow (b, Expr.Num r) -> binomial_of (expand b v n) r n
  | Expr.Pow _ -> raise (Cannot_expand "symbolic exponent")
  | Expr.Fun (f, args) -> fn f (List.map (fun a -> expand a v n) args) n

And the function cases.

and fn (f : string) (args : Series.t list) (n : int) : Series.t =
  let one = Series.const Rational.one n in
  match (f, args) with
  | "exp", [ u ] -> Series.exp_s u
  | "log", [ u ] -> Series.log1p_s (unit_part u n)
  | "sin", [ u ] -> Series.sin_s u
  | "cos", [ u ] -> Series.cos_s u
  | "tan", [ u ] -> Series.div (Series.sin_s u) (Series.cos_s u)
  | "sqrt", [ u ] -> binomial_of u half n
  | "sinh", [ u ] ->
    Series.scale (Series.sub (Series.exp_s u) (Series.exp_s (Series.neg u))) half
  | "cosh", [ u ] ->
    Series.scale (Series.add (Series.exp_s u) (Series.exp_s (Series.neg u))) half
  | "tanh", [ u ] ->
    let e = Series.exp_s (Series.scale u (Rational.of_int 2)) in
    Series.div (Series.sub e one) (Series.add e one)
  (* atan u = integral of u' / (1 + u^2) *)
  | "atan", [ u ] -> Series.integrate (Series.div (Series.diff u) (Series.add one (Series.mul u u)))
  | _ -> raise (Cannot_expand ("no expansion for " ^ f))

§ 07

Limits

Once the expansion exists the limit is not a computation, it is a look at the valuation.

Three cases, and nothing else.

(* Read the limit off the valuation: a positive valuation means every
   term vanishes, zero means the constant term is the answer, and a
   negative one means it blows up with the sign of the leading
   coefficient. *)
let of_series (s : Series.t) : result =
  let s = Series.normalize s in
  if Series.is_zero s then Finite Expr.num_zero
  else if s.Series.v > 0 then Finite Expr.num_zero
  else if s.Series.v = 0 then Finite (Expr.Num (Series.coeff s 0))
  else if Rational.sign (Series.coeff s s.Series.v) > 0 then PosInf
  else NegInf

A leading coefficient that cancels at one truncation may be nonzero at the next, so the order is increased rather than trusted.

(* Try increasing orders: a leading coefficient that cancels to zero at
   one truncation may be nonzero at the next, and only a series that is
   genuinely zero stays zero. *)
let with_orders (f : int -> 'a) : 'a =
  let rec go = function
    | [] -> raise (Cannot_expand "no order sufficed")
    | n :: rest -> ( try f n with Cannot_expand _ when rest <> [] -> go rest)
  in
  go [ 10; 20; 32 ]

x -> infinity is the same machinery under x = 1/t.

(* x -> infinity becomes t -> 0 under x = 1/t. *)
let at_infinity (e : Expr.t) (v : string) : result =
  with_orders (fun n ->
    let sub = Expr.substitute v (Expr.Pow (Expr.Sym v, Expr.int (-1))) (Expr.simplify e) in
    of_series (expand sub v n))

Verbatim.

> limit sin(x)/x, x, 0
1
> limit (1 - cos(x))/x^2, x, 0
1/2
> limit (exp(x) - 1 - x)/x^2, x, 0
1/2
> limit 1/x^2, x, 0
+infinity
> limit (2*x^2 + 3)/(x^2 - 1), x, oo
2
> limit atan(1/x), x, oo
0

This is not Gruntz. There is no most-rapidly-varying subexpression analysis, so a limit whose answer is invisible in any finite truncation, anything needing an exponential tower, is refused rather than guessed at.

§ 08

Integration, and what is decidable

Differentiation is total and integration is not, which changes what the code can promise.

differentiation    total      every elementary function has an
                              elementary derivative
 
integration        partial    exp(x^2) has no elementary antiderivative,
                              and deciding that in general is the Risch
                              algorithm, which is not here
 
rational functions decidable  every one has an elementary antiderivative,
                              and it is a rational part plus logarithms
                              and arctangents. That part is a decision
                              procedure and is implemented in full.

§ 09

Rational functions

Partial fractions from Part 4 reduce the problem to one term at a time. A linear factor gives a logarithm or a power; an irreducible quadratic gives a logarithm and an arctangent.

One term. The numerator is split into a multiple of the denominator's derivative plus a constant, so the first piece integrates by the power rule and only the constant needs the recurrence.

(* One partial-fraction term: numerator / g^j with deg numerator < deg g. *)
let term_integral (ctx : Algebra.ctx) (num : Poly.t) (g : Poly.t) (j : int) (v : string) : Expr.t =
  let x = Expr.Sym v in
  let dg = Poly.degree_in g v in
  let rat p = match Poly.to_rational p with Some r -> r | None -> raise (Cannot_integrate "non-constant coefficient") in
  if dg = 1 then begin
    (* g = a x + b, numerator is a constant k. *)
    let a = rat (Poly.coeff_in g v 1) and b = rat (Poly.coeff_in g v 0) in
    let k = rat num in
    let ge = Expr.Add [ Expr.Mul [ Expr.Num a; x ]; Expr.Num b ] in
    if j = 1 then Expr.simplify (Expr.Mul [ Expr.Num (Rational.div k a); Expr.Fun ("log", [ ge ]) ])
    else
      let c = Rational.div k (Rational.mul a (Rational.of_int (1 - j))) in
      Expr.simplify (Expr.Mul [ Expr.Num c; Expr.Pow (ge, Expr.int (1 - j)) ])
  end
  else if dg = 2 then begin
    let a = rat (Poly.coeff_in g v 2) and b = rat (Poly.coeff_in g v 1) and c = rat (Poly.coeff_in g v 0) in
    let p = rat (Poly.coeff_in num v 1) and q = rat (Poly.coeff_in num v 0) in
    (* Split the numerator into a multiple of g' plus a constant, so the
       first piece integrates by the power rule and only the constant
       needs the recurrence. *)
    let alpha = Rational.div p (Rational.mul (Rational.of_int 2) a) in
    let beta = Rational.sub q (Rational.mul alpha b) in
    let ge = Algebra.of_poly ctx g in
    let first =
      if Rational.is_zero alpha then Expr.num_zero
      else if j = 1 then Expr.Mul [ Expr.Num alpha; Expr.Fun ("log", [ ge ]) ]
      else
        Expr.Mul
          [ Expr.Num (Rational.div alpha (Rational.of_int (1 - j))); Expr.Pow (ge, Expr.int (1 - j)) ]
    in
    let second =
      if Rational.is_zero beta then Expr.num_zero
      else Expr.Mul [ Expr.Num beta; inv_quadratic a b c j v ]
    in
    Expr.simplify (Expr.Add [ first; second ])
  end
  else raise (Cannot_integrate "irreducible factor of degree 3 or more")

The recurrence for the constant part is the standard reduction

\int \frac{d x}{q ^{j}} = \frac{2 a x + b}{( j - 1 ) D q ^{j - 1}} + \frac{2 ( 2 j - 3 ) a}{( j - 1 ) D} \int \frac{d x}{q ^{j - 1}}, D = 4 a c - b^{2}

bottoming out at j = 1, which is the arctangent when D > 0 and a pair of logarithms when it is negative. Both cases occur: x^2 - 2 is irreducible over the rationals and still has real roots.

Which is why the discriminant is tested rather than assumed.

let rec inv_quadratic (a : Rational.t) (b : Rational.t) (c : Rational.t) (j : int) (v : string) : Expr.t =
  let x = Expr.Sym v in
  let q = Expr.Add [ Expr.Mul [ Expr.Num a; Expr.Pow (x, Expr.int 2) ]; Expr.Mul [ Expr.Num b; x ]; Expr.Num c ] in
  let d = Rational.sub (Rational.mul (Rational.of_int 4) (Rational.mul a c)) (Rational.mul b b) in
  if j = 1 then begin
    if Rational.is_zero d then raise (Cannot_integrate "degenerate quadratic");
    if Rational.sign d > 0 then begin
      (* 2/sqrt(D) * atan((2ax + b)/sqrt(D)) *)
      let s = sqrt_expr d in
      Expr.simplify
        (Expr.Mul
           [ Expr.int 2; Expr.Pow (s, Expr.int (-1));
             Expr.Fun ("atan",
               [ Expr.Mul
                   [ Expr.Add [ Expr.Mul [ Expr.Num (Rational.mul (Rational.of_int 2) a); x ]; Expr.Num b ];
                     Expr.Pow (s, Expr.int (-1)) ] ]) ])
    end
    else begin
      (* Real distinct roots: 1/(a(x-r1)(x-r2)) splits into two logs. *)
      let s = sqrt_expr (Rational.neg d) in
      let two_a = Expr.Num (Rational.mul (Rational.of_int 2) a) in
      let u = Expr.Add [ Expr.Mul [ two_a; x ]; Expr.Num b ] in
      Expr.simplify
        (Expr.Mul
           [ Expr.Pow (s, Expr.int (-1));
             Expr.Add
               [ Expr.Fun ("log", [ Expr.Add [ u; Expr.Mul [ Expr.int (-1); s ] ] ]);
                 Expr.Mul [ Expr.int (-1); Expr.Fun ("log", [ Expr.Add [ u; s ] ]) ] ] ])
    end
  end
  else begin
    if Rational.is_zero d then raise (Cannot_integrate "degenerate quadratic");
    let jm = Rational.of_int (j - 1) in
    let denom = Rational.mul jm d in
    let first =
      Expr.Mul
        [ Expr.Add [ Expr.Mul [ Expr.Num (Rational.mul (Rational.of_int 2) a); x ]; Expr.Num b ];
          Expr.Pow (Expr.Num denom, Expr.int (-1));
          Expr.Pow (q, Expr.int (-(j - 1))) ]
    in
    let coeff =
      Rational.div (Rational.mul (Rational.of_int (2 * ((2 * j) - 3))) a) denom
    in
    Expr.simplify (Expr.Add [ first; Expr.Mul [ Expr.Num coeff; inv_quadratic a b c (j - 1) v ] ])
  end

§ 10

The kernel trap

The polynomial layer turns anything non-polynomial into an opaque variable, which is what makes expand work on sin(x)*(x+1)^2. In an integrator it is a bug waiting to happen, and it happened.

What the first version returned. It is not a bad answer, it is a wrong one.

> integrate x*exp(x), x
1/2*exp(x)*x^2          <- wrong
 
exp(x) became an opaque variable, the integrator saw a constant
times x, and integrated it as one. The real answer needs
integration by parts, which is not implemented.

The guard. A kernel that mentions the variable of integration is not a constant, and refusing is the only correct response.

(* The polynomial layer turns anything non-polynomial into an opaque
   variable. That is exactly wrong here if the kernel depends on the
   variable of integration: exp(x) would be carried along as a constant
   and x*exp(x) would integrate to x^2 exp(x)/2. Refuse instead. *)
let check_kernels (ctx : Algebra.ctx) (p : Poly.t) (v : string) : unit =
  List.iter
    (fun name ->
      if name <> v then
        match List.assoc_opt name ctx.Algebra.kernels with
        | Some k when List.mem v (Expr.free_symbols k) ->
          raise (Cannot_integrate ("cannot integrate in terms of " ^ Expr.to_string k))
        | _ -> ())
    (Poly.vars_of p)

§ 11

Verification

Every integral here is checked by differentiating it. That is a real test rather than a restatement: differentiation lives in another module, uses another method, and is decidable where integration is not.

The property.

(* The check that makes the whole thing testable: differentiate the
   answer and compare with the integrand. *)
let verify (e : Expr.t) (v : string) : bool =
  match integrate e v with
  | exception Cannot_integrate _ -> false
  | anti ->
    let d = Deriv.diff anti v in
    Expr.compare_expr (Algebra.together (Expr.Add [ d; Expr.Mul [ Expr.int (-1); e ] ])) Expr.num_zero = 0

Every rational function in the suite, verified this way. Verbatim.

integrand                 antiderivative                                  check
x^2                       1/3*x^3                                         verified
1/x                       log(x)                                          verified
1/(x^2 - 1)               -1/2*log(1 + x) + 1/2*log(-1 + x)               verified
1/(x^2 + 1)               atan(x)                                         verified
1/(x^3 + x)               log(x) - 1/2*log(1 + x^2)                       verified
1/((x - 1)^2*(x + 2))     -1/3*(-1 + x)^(-1) - 1/9*log(-1 + x)
                            + 1/9*log(2 + x)                              verified
1/(x^2 + x + 1)           2*atan(3^(-1/2)*(1 + 2*x))*3^(-1/2)             verified
1/(x^2 + 1)^2             1/2*x*(1 + x^2)^(-1) + 1/2*atan(x)              verified
exp(2*x + 1)              1/2*exp(1 + 2*x)                                verified
log(x)                    -x + x*log(x)                                   verified

§ 12

What is refused, and why that is the feature

Three refusals, all honest.

> integrate x*exp(x), x
cannot integrate: cannot integrate in terms of exp(x)
> integrate sin(x^2), x
cannot integrate: function of a non-linear argument
> integrate tan(x), x
cannot integrate: no table entry for tan

The last one is deliberate and is worth explaining. -log(cos x) is the right answer, and it was in the table until it failed its own verification: nothing in this system knows that tan = sin/cos, so the check could not confirm it. An entry whose correctness the system cannot demonstrate does not belong in the table.

§ 13

Tests

The randomized one: build a rational function with a planted denominator, integrate it, differentiate the result, and demand zero.

(* Randomized: build a rational function with a planted denominator and
   demand that it integrates and differentiates back. *)
let test_random_rational () =
  Random.init 5150;
  let x = Expr.Sym "x" in
  for _ = 1 to 60 do
    let lin () = Expr.Add [ x; Expr.int (Random.int 7 - 3) ] in
    let quad () = Expr.Add [ Expr.Pow (x, Expr.int 2); Expr.int (1 + Random.int 4) ] in
    let den =
      match Random.int 3 with
      | 0 -> Expr.Mul [ lin (); lin () ]
      | 1 -> Expr.Mul [ lin (); quad () ]
      | _ -> quad ()
    in
    let num = Expr.Add [ Expr.Mul [ Expr.int (1 + Random.int 4); x ]; Expr.int (Random.int 5) ] in
    let f = Expr.simplify (Expr.Mul [ num; Expr.Pow (den, Expr.int (-1)) ]) in
    match Integrate.integrate f "x" with
    | exception Integrate.Cannot_integrate _ -> ()
    | exception Division_by_zero -> ()
    | anti ->
      let back = Algebra.together (Expr.Add [ Deriv.diff anti "x"; Expr.Mul [ Expr.int (-1); f ] ]) in
      check "a random rational function integrates back"
        (Expr.compare_expr back Expr.num_zero = 0)
  done

§ 14

What is deliberately not here

Three gaps, in order of how much they hurt.

Risch                 the decision procedure for elementary integration.
                      Without it, integrate refuses rather than proving
                      that exp(x^2) has no elementary antiderivative.
 
integration by parts  and substitution. x*exp(x) is refused, not solved.
 
Gruntz                limits are read off a series, so anything needing
                      an exponential tower is out of reach.

§ 15

Download

The snapshot at the end of this part. New in Part 5: deriv.ml, series.ml, limit.ml, integrate.ml, and the tests test_deriv.ml, test_series.ml, test_limit.ml, test_integrate.ml. Unchanged from Part 4: bigint.ml, rational.ml, expr.ml, fp.ml, upoly.ml, poly.ml, factor.ml, ratfun.ml, parser.ml, and main.ml.

← older

Building a CAS in OCaml Part 4

2026-04-27

newer →

async in Haskell

2026-05-24

Building a CAS in OCaml Part 5

What Part 3 refused

The general power rule

Laurent series

Elementary functions, by recurrence

Where the known order goes

Expanding an expression

Limits

Integration, and what is decidable

Rational functions

The kernel trap

Verification

What is refused, and why that is the feature

Tests

What is deliberately not here

Download

Building a CAS in OCaml Part 4

Building a CAS in OCaml Part 3

Building a CAS in OCaml Part 2

Building a CAS in OCaml Part 4

async in Haskell