Building a CAS in OCaml Part 1

2026-04-06 · 13 min

A computer algebra system is exact or it is nothing, and exactness bottoms out in integers that do not overflow. This is the whole bignum layer: representation, addition, two multiplication algorithms, long division, and decimal conversion - in OCaml, with no dependencies, tested against native arithmetic where it fits and against algebraic identities where it does not.

§ 01

Layout, and how to build it

The snapshot this post produces. Every later part of the series adds to the same lib/ and keeps its own complete copy, so a diff between two parts is exactly what that post changed.

ocaml-cas/part1/
├── dune-project
├── build.sh
├── lib/
│   ├── dune
│   └── bigint.ml
├── bin/
│   ├── dune
│   └── main.ml
└── test/
    ├── dune
    └── test_bigint.ml

dune-project, lib/dune, bin/dune, test/dune - the library is named cas, so the module below is Cas.Bigint.

(lang dune 3.16)
 
(library
 (name cas))
 
(executable
 (name main)
 (libraries cas))
 
(test
 (name test_bigint)
 (libraries cas))

§ 02

Representation

Sign and magnitude, with the magnitude a little-endian array of limbs in base 2^30 and no leading zero limbs. One representation per value, so structural comparison of magnitudes is meaningful.

Base 2^30 is the interesting choice: a product of two limbs is under 2^60, and an OCaml native int holds up to 2^62 - 1, which leaves room for a carry and an addend without ever reaching for Int64 or boxing.

let base_bits = 30
let base = 1 lsl base_bits
let base_mask = base - 1

Strip leading zero limbs so every magnitude is canonical. The identity check avoids a copy in the common case where nothing needed stripping.

let mag_normalize (a : int array) : int array =
  let n = ref (Array.length a) in
  while !n > 0 && a.(!n - 1) = 0 do
    decr n
  done;
  if !n = Array.length a then a else Array.sub a 0 !n
 
let mag_is_zero (a : int array) : bool = Array.length a = 0

Because magnitudes are canonical, a longer array is always the larger number - only equal lengths need a limb-by-limb walk, from the top down.

let mag_compare (a : int array) (b : int array) : int =
  let la = Array.length a and lb = Array.length b in
  if la <> lb then compare la lb
  else begin
    let rec go i =
      if i < 0 then 0
      else if a.(i) <> b.(i) then compare a.(i) b.(i)
      else go (i - 1)
    in
    go (la - 1)
  end

§ 03

Addition and subtraction of magnitudes

Carry propagation, one limb at a time, with the top limb of the result holding the final carry before normalization drops it if it is zero.

let mag_add (a : int array) (b : int array) : int array =
  let la = Array.length a and lb = Array.length b in
  let l = max la lb in
  let r = Array.make (l + 1) 0 in
  let carry = ref 0 in
  for i = 0 to l - 1 do
    let s = (if i < la then a.(i) else 0) + (if i < lb then b.(i) else 0) + !carry in
    r.(i) <- s land base_mask;
    carry := s lsr base_bits
  done;
  r.(l) <- !carry;
  mag_normalize r

Subtraction assumes a >= b, which the signed layer guarantees by comparing first and swapping. Borrowing is the mirror of carrying.

(* Requires mag_compare a b >= 0. *)
let mag_sub (a : int array) (b : int array) : int array =
  let la = Array.length a and lb = Array.length b in
  let r = Array.make la 0 in
  let borrow = ref 0 in
  for i = 0 to la - 1 do
    let d = a.(i) - (if i < lb then b.(i) else 0) - !borrow in
    if d < 0 then begin
      r.(i) <- d + base;
      borrow := 1
    end
    else begin
      r.(i) <- d;
      borrow := 0
    end
  done;
  mag_normalize r

Shifting by whole limbs is multiplication by a power of the base - needed by Karatsuba. Shifting by bits within a limb is needed by division.

let mag_shift_limbs (a : int array) (k : int) : int array =
  if mag_is_zero a || k = 0 then a
  else begin
    let n = Array.length a in
    let r = Array.make (n + k) 0 in
    Array.blit a 0 r k n;
    r
  end
 
(* 0 <= s < base_bits *)
let mag_shift_left_bits (a : int array) (s : int) : int array =
  if s = 0 || mag_is_zero a then a
  else begin
    let n = Array.length a in
    let r = Array.make (n + 1) 0 in
    let carry = ref 0 in
    for i = 0 to n - 1 do
      let t = (a.(i) lsl s) lor !carry in
      r.(i) <- t land base_mask;
      carry := t lsr base_bits
    done;
    r.(n) <- !carry;
    mag_normalize r
  end
 
(* 0 <= s < base_bits *)
let mag_shift_right_bits (a : int array) (s : int) : int array =
  if s = 0 || mag_is_zero a then a
  else begin
    let n = Array.length a in
    let r = Array.make n 0 in
    let carry = ref 0 in
    for i = n - 1 downto 0 do
      let cur = a.(i) in
      r.(i) <- (cur lsr s) lor (!carry lsl (base_bits - s));
      carry := cur land ((1 lsl s) - 1)
    done;
    mag_normalize r
  end

§ 04

Schoolbook multiplication

The classic double loop, accumulating into the result in place. Skipping a zero multiplier limb is worth the branch on sparse-ish inputs.

let mag_mul_school (a : int array) (b : int array) : int array =
  let la = Array.length a and lb = Array.length b in
  if la = 0 || lb = 0 then [||]
  else begin
    let r = Array.make (la + lb) 0 in
    for i = 0 to la - 1 do
      let ai = a.(i) in
      if ai <> 0 then begin
        let carry = ref 0 in
        for j = 0 to lb - 1 do
          (* ai * b.(j) < 2^60, plus two addends < 2^31: no overflow. *)
          let t = (ai * b.(j)) + r.(i + j) + !carry in
          r.(i + j) <- t land base_mask;
          carry := t lsr base_bits
        done;
        let k = ref (i + lb) in
        while !carry <> 0 do
          let t = r.(!k) + !carry in
          r.(!k) <- t land base_mask;
          carry := t lsr base_bits;
          incr k
        done
      end
    done;
    mag_normalize r
  end

The overflow argument, written out, because it is the reason base 2^30 was picked and the thing that silently breaks if someone raises it to 2^31.

ai        <= 2^30 - 1
b.(j)     <= 2^30 - 1
ai*b.(j)  <= 2^60 - 2^31 + 1
r.(i+j)   <= 2^30 - 1
carry     <= 2^30
sum       <  2^60 + 2^31  <  2^62 - 1 = max_int
 
at base_bits = 31 the product alone reaches 2^62 and this overflows

§ 05

Karatsuba

Splitting each operand in half turns one big multiplication into four half-sized ones; the trick is that the middle term can be recovered from a single extra product rather than two, taking it from four to three and the exponent from log2(4) = 2 down to log2(3) ~ 1.585.

The identity, before the code.

a = a1*B^m + a0
b = b1*B^m + b0
 
a*b = a1*b1 * B^2m
    + (a1*b0 + a0*b1) * B^m
    + a0*b0
 
and the middle term is recoverable without computing it directly:
  a1*b0 + a0*b1 = (a0 + a1)*(b0 + b1) - a0*b0 - a1*b1
 
so three products (a0*b0, a1*b1, (a0+a1)*(b0+b1)) suffice, not four

Recursion bottoming out in the schoolbook version, which wins below roughly thirty limbs where the extra additions and allocations cost more than the saved product.

let karatsuba_threshold = 32
 
let rec mag_mul (a : int array) (b : int array) : int array =
  let la = Array.length a and lb = Array.length b in
  if la = 0 || lb = 0 then [||]
  else if la < karatsuba_threshold || lb < karatsuba_threshold then
    mag_mul_school a b
  else begin
    let m = (max la lb + 1) / 2 in
    let low x =
      let n = Array.length x in
      if n <= m then x else mag_normalize (Array.sub x 0 m)
    in
    let high x =
      let n = Array.length x in
      if n <= m then [||] else mag_normalize (Array.sub x m (n - m))
    in
    let a0 = low a and a1 = high a in
    let b0 = low b and b1 = high b in
    let z0 = mag_mul a0 b0 in
    let z2 = mag_mul a1 b1 in
    let z1 = mag_sub (mag_sub (mag_mul (mag_add a0 a1) (mag_add b0 b1)) z0) z2 in
    mag_add (mag_add z0 (mag_shift_limbs z1 m)) (mag_shift_limbs z2 (2 * m))
  end

Both mag_sub calls are safe because z1 is mathematically non-negative: (a0+a1)(b0+b1) always dominates a0*b0 + a1*b1 for non-negative limbs. That is a real precondition, not an accident, and the test suite checks it by comparing against the schoolbook result on random inputs.

$ ./_out/main bench 8000

Measured on this machine, random operands, both algorithms on the same inputs with the result compared for equality every time.

limbs    schoolbook   karatsuba   agree
  500       0.000s      0.001s     true
 1000       0.002s      0.002s     true
 2000       0.006s      0.005s     true
 4000       0.024s      0.013s     true
 8000       0.098s      0.035s     true
16000       0.392s      0.105s     true

§ 06

Division: the short case

Dividing by a single limb is a straight walk from the top down, carrying the remainder. The intermediate stays under 2^60 because the running remainder is always below the divisor, which is below the base.

(* 0 < d < base *)
let mag_divmod_small (a : int array) (d : int) : int array * int =
  let n = Array.length a in
  let q = Array.make n 0 in
  let r = ref 0 in
  for i = n - 1 downto 0 do
    let cur = (!r * base) + a.(i) in
    q.(i) <- cur / d;
    r := cur mod d
  done;
  (mag_normalize q, !r)

§ 07

Division: Knuth algorithm D

The general case is the one piece of a bignum library that is genuinely hard to get right. Normalize so the divisor's top limb has its high bit set, estimate each quotient limb from the top two limbs of the running remainder, correct the estimate down by at most two, then multiply-and-subtract - with an add-back path for the rare case where the estimate was still one too large.

Normalization is what bounds the error in the estimate: with the divisor's top limb at least base/2, the two-limb estimate is never more than two above the true quotient limb.

let bit_length (x : int) : int =
  let rec go x acc = if x = 0 then acc else go (x lsr 1) (acc + 1) in
  go x 0
 
let mag_divmod (u : int array) (v : int array) : int array * int array =
  let n = Array.length v in
  if n = 0 then raise Division_by_zero;
  if mag_compare u v < 0 then ([||], u)
  else if n = 1 then begin
    let q, r = mag_divmod_small u v.(0) in
    (q, if r = 0 then [||] else [| r |])
  end
  else begin
    let s = base_bits - bit_length v.(n - 1) in
    let vn = mag_shift_left_bits v s in
    let un0 = mag_shift_left_bits u s in
    let vl = Array.length vn in
    let m = Array.length un0 - vl in
    (* One extra high limb so index j + vl is always in range. *)
    let un = Array.append un0 [| 0 |] in
    let q = Array.make (m + 1) 0 in

The estimate, and its correction loop. qhat starts as a two-limb-by-one-limb quotient and is walked down until the next limb of the divisor is consistent with it.

    for j = m downto 0 do
      let num = (un.(j + vl) * base) + un.(j + vl - 1) in
      let qhat = ref (num / vn.(vl - 1)) in
      let rhat = ref (num mod vn.(vl - 1)) in
      let continue_ = ref true in
      while !continue_ do
        if !qhat >= base || !qhat * vn.(vl - 2) > (!rhat * base) + un.(j + vl - 2)
        then begin
          decr qhat;
          rhat := !rhat + vn.(vl - 1);
          if !rhat >= base then continue_ := false
        end
        else continue_ := false
      done;

Multiply the divisor by the estimate and subtract it from the running remainder in one pass, tracking a multiplication carry and a subtraction borrow simultaneously.

      let borrow = ref 0 and carry = ref 0 in
      for i = 0 to vl - 1 do
        let p = (!qhat * vn.(i)) + !carry in
        carry := p lsr base_bits;
        let t = un.(i + j) - (p land base_mask) - !borrow in
        if t < 0 then begin
          un.(i + j) <- t + base;
          borrow := 1
        end
        else begin
          un.(i + j) <- t;
          borrow := 0
        end
      done;
      let t = un.(j + vl) - !carry - !borrow in
      if t < 0 then begin
        un.(j + vl) <- t + base;
        borrow := 1
      end
      else begin
        un.(j + vl) <- t;
        borrow := 0
      end;

The add-back path. A borrow out of the top limb means the estimate was one too large after all - give one back by adding the divisor straight back in. This branch is famously rare (roughly two in 2^30 random inputs) which is exactly why it is worth having a randomized test that reaches it.

      if !borrow <> 0 then begin
        decr qhat;
        let c = ref 0 in
        for i = 0 to vl - 1 do
          let t = un.(i + j) + vn.(i) + !c in
          un.(i + j) <- t land base_mask;
          c := t lsr base_bits
        done;
        un.(j + vl) <- (un.(j + vl) + !c) land base_mask
      end;
      q.(j) <- !qhat
    done;
    let rem = mag_normalize (Array.sub un 0 vl) in
    (mag_normalize q, mag_shift_right_bits rem s)
  end

§ 08

The signed layer

Sign as an int rather than a variant, so the sign of a product is a multiplication rather than a match. Zero is the only value with sign 0, and make is the only constructor that can produce it.

type t = {
  sign : int;        (* 1, -1, or 0 when the value is zero *)
  mag : int array;   (* little-endian base-2^30 limbs, no leading zeros *)
}
 
let zero = { sign = 0; mag = [||] }
let one = { sign = 1; mag = [| 1 |] }
let minus_one = { sign = -1; mag = [| 1 |] }
 
let make (sign : int) (mag : int array) : t =
  let mag = mag_normalize mag in
  if mag_is_zero mag then zero else { sign; mag }
 
let neg (a : t) : t = if a.sign = 0 then zero else { a with sign = -a.sign }
let abs (a : t) : t = if a.sign = 0 then zero else { a with sign = 1 }

of_int, including the one value that needs care: min_int has no positive counterpart, so negating it to get a magnitude would wrap. A logical shift reads the word as unsigned and sidesteps it entirely.

let of_int (n : int) : t =
  if n = 0 then zero
  else if n > 0 then { sign = 1; mag = mag_of_pos_int n }
  else if n <> min_int then { sign = -1; mag = mag_of_pos_int (-n) }
  else
    {
      sign = -1;
      mag =
        mag_normalize
          [|
            n land base_mask;
            (n lsr base_bits) land base_mask;
            (n lsr (2 * base_bits)) land base_mask;
          |];
    }

Addition dispatches on whether the signs agree; when they differ it is a subtraction of the smaller magnitude from the larger, with the sign of the larger.

let add (a : t) (b : t) : t =
  if a.sign = 0 then b
  else if b.sign = 0 then a
  else if a.sign = b.sign then { sign = a.sign; mag = mag_add a.mag b.mag }
  else begin
    let c = mag_compare a.mag b.mag in
    if c = 0 then zero
    else if c > 0 then make a.sign (mag_sub a.mag b.mag)
    else make b.sign (mag_sub b.mag a.mag)
  end
 
let sub (a : t) (b : t) : t = add a (neg b)
 
let mul (a : t) (b : t) : t =
  if a.sign = 0 || b.sign = 0 then zero
  else make (a.sign * b.sign) (mag_mul a.mag b.mag)

Division truncates toward zero and the remainder takes the sign of the dividend, matching OCaml's own (/) and (mod) rather than a mathematical floor division - the test suite pins this by checking against native ints on every small pair.

let divmod (a : t) (b : t) : t * t =
  if b.sign = 0 then raise Division_by_zero;
  if a.sign = 0 then (zero, zero)
  else begin
    let q, r = mag_divmod a.mag b.mag in
    (make (a.sign * b.sign) q, make a.sign r)
  end
 
let div (a : t) (b : t) : t = fst (divmod a b)
let rem (a : t) (b : t) : t = snd (divmod a b)
 
let rec gcd (a : t) (b : t) : t =
  if b.sign = 0 then abs a else gcd b (rem a b)
 
let pow (a : t) (e : int) : t =
  if e < 0 then invalid_arg "Bigint.pow: negative exponent";
  let rec go acc b e =
    if e = 0 then acc
    else if e land 1 = 1 then go (mul acc b) (mul b b) (e lsr 1)
    else go acc (mul b b) (e lsr 1)
  in
  go one a e

§ 09

Decimal conversion

Printing divides out nine decimal digits at a time, since 10^9 is the largest power of ten a single limb division handles. Every chunk after the first is zero-padded, which is where a naive implementation drops digits.

let decimal_chunk = 1_000_000_000
let decimal_chunk_digits = 9
 
let to_string (a : t) : string =
  if a.sign = 0 then "0"
  else begin
    let rec go m acc =
      if mag_is_zero m then acc
      else
        let q, r = mag_divmod_small m decimal_chunk in
        go q (r :: acc)
    in
    match go a.mag [] with
    | [] -> "0"
    | first :: rest ->
      let buf = Buffer.create 32 in
      if a.sign < 0 then Buffer.add_char buf '-';
      Buffer.add_string buf (string_of_int first);
      List.iter
        (fun c ->
          Buffer.add_string buf
            (Printf.sprintf "%0*d" decimal_chunk_digits c))
        rest;
      Buffer.contents buf
  end

Parsing absorbs nine digits at a time for the same reason, so the number of bignum multiplications is a ninth of the digit count. The first chunk is deliberately the short one, so every subsequent chunk is exactly nine digits wide.

let of_string (s : string) : t =
  let s = String.trim s in
  let len = String.length s in
  if len = 0 then invalid_arg "Bigint.of_string: empty string";
  let negative, start =
    match s.[0] with
    | '-' -> (true, 1)
    | '+' -> (false, 1)
    | _ -> (false, 0)
  in
  if start >= len then invalid_arg "Bigint.of_string: no digits";
  let acc = ref zero in
  let chunk_mul = of_int decimal_chunk in
  let i = ref start in
  let first_len =
    let total = len - start in
    let r = total mod decimal_chunk_digits in
    if r = 0 then decimal_chunk_digits else r
  in
  let read_chunk n =
    let v = ref 0 in
    for _ = 1 to n do
      let c = s.[!i] in
      if c < '0' || c > '9' then
        invalid_arg (Printf.sprintf "Bigint.of_string: bad character %C" c);
      v := (!v * 10) + (Char.code c - Char.code '0');
      incr i
    done;
    !v
  in
  acc := of_int (read_chunk first_len);
  while !i < len do
    let chunk = read_chunk decimal_chunk_digits in
    acc := add (mul !acc chunk_mul) (of_int chunk)
  done;
  if negative then neg !acc else !acc

§ 10

Testing something that cannot be checked by eye

Two strategies, because neither alone is enough. Where values fit in a native int, check against native arithmetic directly. Where they do not, check algebraic identities that must hold regardless of magnitude.

Against native ints, exhaustively over a grid of interesting values - including the limb boundary at 2^30 and both extremes of the native range.

let test_divmod_small () =
  let vals = [ 1; -1; 3; -3; 7; -7; 1000; -1000; 1073741823; -1073741824 ] in
  let nums = [ 0; 1; -1; 20; -20; 999999999; -999999999; max_int ] in
  List.iter
    (fun a ->
      List.iter
        (fun b ->
          let q, r = divmod (of_int a) (of_int b) in
          check_eq (Printf.sprintf "%d / %d" a b) (of_int (a / b)) q;
          check_eq (Printf.sprintf "%d mod %d" a b) (of_int (a mod b)) r)
        vals)
    nums

Beyond native range, identities: a = q*b + r with |r| < |b| is the one that actually catches algorithm D bugs, including the add-back path.

let test_random_division () =
  Random.init 424242;
  for _ = 1 to 200 do
    let la = 1 + Random.int 40 and lb = 1 + Random.int 20 in
    let a = make 1 (random_mag la) in
    let b = make 1 (random_mag lb) in
    if not (is_zero b) then begin
      let q, r = divmod a b in
      check_eq "random division identity" a (add (mul q b) r);
      check "random remainder in range"
        (compare r zero >= 0 && compare r b < 0)
    end
  done

Karatsuba is checked against the algorithm it replaces rather than against expected values - the strongest available oracle, since the schoolbook version is simple enough to trust.

(* Random.bits () yields exactly 30 random bits, which is exactly one
   limb in base 2^30. (Random.int base would raise: its bound must be
   strictly below 2^30.) *)
let random_mag n = Array.init n (fun _ -> Random.bits ())
 
let test_karatsuba_agrees () =
  Random.init 20260919;
  for _ = 1 to 50 do
    let la = 1 + Random.int 120 and lb = 1 + Random.int 120 in
    let a = mag_normalize (random_mag la) in
    let b = mag_normalize (random_mag lb) in
    let school = mag_mul_school a b in
    let karatsuba = mag_mul a b in
    incr checks;
    if mag_compare school karatsuba <> 0 then begin
      incr failures;
      Printf.printf "FAIL: karatsuba disagrees on %d x %d limbs\n" la lb
    end
  done

A known-answer test worth having, because it exercises hundreds of multiplications and a long decimal conversion at once: 100! has 158 digits, begins 933262154439, and ends in exactly 24 zeros.

  let rec fact n acc = if n <= 1 then acc else fact (n - 1) (mul acc (of_int n)) in
  let f100 = to_string (fact 100 one) in
  check "100! digit count" (String.length f100 = 158);
  check "100! prefix"
    (String.length f100 >= 12 && String.sub f100 0 12 = "933262154439");
  check "100! trailing zeros" (trailing_zeros f100 = 24)

§ 11

The driver

Enough of a command line to poke at the library by hand.

$ ./_out/main fact 30
265252859812191058636308480000000
$ ./_out/main pow 2 100
1267650600228229401496703205376
$ ./_out/main gcd 123456789012345678901234567890 987654321098765432109876543210
9000000000900000000090
$ ./_out/main div 100000000000000000000000000000000 7
quotient  14285714285714285714285714285714
remainder 2

§ 12

Download

The complete snapshot as it stands at the end of this part: bigint.ml, test_bigint.ml, main.ml. Part 2 builds rationals and the expression tree on exactly this module, and keeps its own full copy so the diff between parts is the change that post describes.

referenced by

Building a CAS in OCaml Part 2
2026-04-13

newer →

Building a CAS in OCaml Part 2

2026-04-13

Building a CAS in OCaml Part 1

Layout, and how to build it

Representation

Addition and subtraction of magnitudes

Schoolbook multiplication

Karatsuba

Division: the short case

Division: Knuth algorithm D

The signed layer

Decimal conversion

Testing something that cannot be checked by eye

The driver

Download

Building a CAS in OCaml Part 2

Building a CAS in OCaml Part 5

Building a CAS in OCaml Part 4

Building a CAS in OCaml Part 3

Building a CAS in OCaml Part 2