Knuth's Algorithm D

Knuth's Algorithm D divides one multiple-precision integer by another. It is the long division taught in school, carried out in base $β = 2^{k}$ on machine words called limbs, with one difficulty that the school version hides: each digit of the quotient has to be guessed. Algorithm D guesses it from the two leading limbs of the current remainder and the leading limb of the divisor. After the divisor is normalized so that its leading limb is at least $β /2$ , the guess is never too small and at most two too large; a test involving one more limb removes almost every error, and a rare add-back step corrects what remains^[1]. Dividing an $(m + n)$ -limb number by an $n$ -limb number takes $O (mn)$ limb operations. The algorithm is given in the second volume of Knuth's The Art of Computer Programming, and is the division used, with variations, by most arbitrary-precision arithmetic libraries.

§ 01

Overview

A multiple-precision integer is stored as a sequence of limbs, $u = \sum_{i} u_{i} β^{i}$ with $0 \leq u_{i} < β$ . Addition, subtraction and multiplication by the schoolbook methods need only single-limb operations and carries. Division is harder because long division needs, at each step, the largest digit $q_{j}$ such that $q_{j} \cdot v$ does not exceed the current partial remainder. In base 10 a person finds it by inspection and trial; in base $2^{28}$ or $2^{64}$ trial and error over $β$ candidates is out of the question. Algorithm D computes a close estimate with one double-limb division and corrects it in at most a few steps.

§ 02

The algorithm

Let $u = (u_{m + n - 1} \dots u_{1} u_{0})_{β}$ and $v = (v_{n - 1} \dots v_{1} v_{0})_{β}$ with $v_{n - 1} \neq = 0$ and $n \geq 2$ . The algorithm computes $q = ⌊ u / v ⌋$ and $r = u mod v$ in these steps^[1]:

D1. Normalize: multiply $u$ and $v$ by $d = 2^{s}$ , chosen so that the leading limb of $v$ satisfies $v_{n - 1} \geq β /2$ . $u$ gains a leading limb $u_{m + n}$ , possibly zero. The quotient is unchanged.

D2–D7. For $j = m, m - 1, \dots, 0$ , divide the $(n + 1)$ -limb number $(u_{j + n} \dots u_{j})_{β}$ by $v$ to obtain the quotient digit $q_{j}$ , and replace those limbs of $u$ by the remainder: estimate $q^$ , refine it, multiply $v$ by $q^$ and subtract, and if the result is negative decrease $q^$ by one and add $v$ back.

D8. Unnormalize: the remainder is the low $n$ limbs of $u$ divided by $d$ .

§ 03

Estimating the quotient digit

The estimate uses the two leading limbs of the partial remainder and the leading limb of the divisor:

q^= min (⌊ \frac{u _{j + n} β + u _{j + n - 1}}{v _{n - 1}} ⌋, β - 1)

Two facts make this work^[1]. The estimate is never too small, $q^\geq q_{j}$ , because the partial remainder is less than $β v$ and the ignored limbs can only make the true quotient smaller. And if $v_{n - 1} \geq ⌊ β /2 ⌋$ , the estimate is at most two too large:

v_{n - 1} \geq ⌊ \frac{β}{2} ⌋ ⟹ q^- 2 \leq q_{j} \leq q^

Without normalization the error is unbounded in terms of the base: with $v_{n - 1} = 1$ the estimate ignores almost all of the divisor. Normalizing by a power of two is a shift and does not change the quotient, which is why it is done first. Pope and Stein had earlier analysed the effect of the leading digit of the divisor on such estimates^[2].

Refinement

A second test with the next limbs of both numbers removes most of the remaining error before any multiplication is done. With $r^= (u_{j + n} β + u_{j + n - 1}) - q^v_{n - 1}$ :

while q^\geq β or q^v_{n - 2} > r^β + u_{j + n - 2} : q^\leftarrow q^- 1, r^\leftarrow r^+ v_{n - 1}

stopping as soon as $r^\geq β$ , since the test can then no longer succeed. It amounts to estimating with three limbs of the remainder and two of the divisor, and it removes every case in which $q^$ is two too large and most in which it is one too large^[1]. What is left is detected in the next step: if subtracting $q^v$ makes the partial remainder negative, $q^$ was one too large, and adding $v$ back restores it.

§ 04

Implementation

Natural numbers as arrays of limbs, with schoolbook addition and multiplication.

(* Natural numbers as little-endian arrays of limbs in base 2^bits, with no
   leading zero limbs. Everything except division is schoolbook. *)
module Make (B : sig val bits : int end) = struct
  let bits = B.bits
  let base = 1 lsl bits
  let mask = base - 1
 
  let trim a = let n = ref (Array.length a) in while !n > 0 && a.(!n - 1) = 0 do decr n done; Array.sub a 0 !n
  let of_int x = let rec go x acc = if x = 0 then acc else go (x lsr bits) ((x land mask) :: acc) in Array.of_list (List.rev (go x []))
 
  let compare a b =
    let la = Array.length a and lb = Array.length b in
    if la <> lb then Stdlib.compare la lb
    else let rec go i = if i < 0 then 0 else if a.(i) <> b.(i) then Stdlib.compare a.(i) b.(i) else go (i - 1) in go (la - 1)
 
  let add a b =
    let n = max (Array.length a) (Array.length b) in
    let r = Array.make (n + 1) 0 and c = ref 0 in
    for i = 0 to n - 1 do
      let s = (if i < Array.length a then a.(i) else 0) + (if i < Array.length b then b.(i) else 0) + !c in
      r.(i) <- s land mask; c := s lsr bits
    done;
    r.(n) <- !c; trim r
 
  let mul a b =
    let r = Array.make (Array.length a + Array.length b + 1) 0 in
    Array.iteri (fun i x ->
      let c = ref 0 in
      Array.iteri (fun j y -> let t = r.(i + j) + x * y + !c in r.(i + j) <- t land mask; c := t lsr bits) b;
      let k = ref (i + Array.length b) in
      while !c > 0 do let t = r.(!k) + !c in r.(!k) <- t land mask; c := t lsr bits; incr k done) a;
    trim r
 
  (* Multiply by a single limb and add a single limb: used for parsing. *)
  let mul_small_add a m d = add (mul a (of_int m)) (of_int d)
  let of_string s = String.fold_left (fun acc ch -> mul_small_add acc 10 (Char.code ch - 48)) [||] s
end

Algorithm D, with counters for how far each first estimate was from the true digit and for add-backs.

(* Knuth's Algorithm D: divide u (m + n limbs) by v (n >= 2 limbs). *)
module Make (B : sig val bits : int end) = struct
  include Nat.Make (B)
 
  type stats = { mutable qhat_over : int array; mutable add_backs : int; mutable digits : int }
  let stats = { qhat_over = Array.make 3 0; add_backs = 0; digits = 0 }
 
  let shift_left a s = (* by s < bits bits, keeping one extra limb *)
    let n = Array.length a in
    let r = Array.make (n + 1) 0 in
    for i = n - 1 downto 0 do
      r.(i + 1) <- r.(i + 1) lor (a.(i) lsr (bits - s));
      r.(i) <- (a.(i) lsl s) land mask
    done;
    if s = 0 then r.(n) <- 0;
    r
 
  let divmod_small u d = (* one-limb divisor *)
    let q = Array.make (Array.length u) 0 and r = ref 0 in
    for i = Array.length u - 1 downto 0 do
      let t = (!r lsl bits) lor u.(i) in q.(i) <- t / d; r := t mod d
    done;
    (trim q, of_int !r)
 
  let divmod u v =
    let n = Array.length v in
    if n = 0 then raise Division_by_zero
    else if compare u v < 0 then ([||], u)
    else if n = 1 then divmod_small u v.(0)
    else begin
      (* D1. Normalize: shift so the divisor's top limb is at least base/2. *)
      let s = let rec go s x = if x >= base / 2 then s else go (s + 1) (x lsl 1) in go 0 v.(n - 1) in
      let v = Array.sub (shift_left v s) 0 n and u = shift_left u s in
      let m = Array.length u - n - 1 in
      let q = Array.make (m + 1) 0 in
      for j = m downto 0 do
        (* D3. Estimate the quotient digit from the top two limbs of the
           current remainder and the top limb of the divisor. *)
        let num = (u.(j + n) lsl bits) + u.(j + n - 1) in
        let qhat = ref (num / v.(n - 1)) and rhat = ref (num mod v.(n - 1)) in
        (* Record how far the estimate is from the true digit, for the
           statistics; the true digit is computed at the end. *)
        let estimate = min !qhat (base - 1) in
        (* Correct with the next limb: at most two decrements. *)
        let testing = ref true in
        while !testing && (!qhat >= base || !qhat * v.(n - 2) > (!rhat lsl bits) + u.(j + n - 2)) do
          decr qhat;
          rhat := !rhat + v.(n - 1);
          if !rhat >= base then testing := false (* the test can no longer succeed *)
        done;
        (* D4. Multiply and subtract qhat * v from u[j .. j+n]. *)
        let borrow = ref 0 and carry = ref 0 in
        for i = 0 to n do
          let p = (if i < n then !qhat * v.(i) else 0) + !carry in
          carry := p lsr bits;
          let t = u.(i + j) - (p land mask) - !borrow in
          if t < 0 then (u.(i + j) <- t + base; borrow := 1) else (u.(i + j) <- t; borrow := 0)
        done;
        (* D6. If the result went negative, qhat was one too large: add v back. *)
        if !borrow <> 0 then begin
          stats.add_backs <- stats.add_backs + 1;
          decr qhat;
          let c = ref 0 in
          for i = 0 to n do
            let t = u.(i + j) + (if i < n then v.(i) else 0) + !c in
            u.(i + j) <- t land mask; c := t lsr bits
          done
        end;
        q.(j) <- !qhat;
        stats.digits <- stats.digits + 1;
        let over = estimate - !qhat in
        if over <= 2 then stats.qhat_over.(over) <- stats.qhat_over.(over) + 1
      done;
      (* D8. Unnormalize the remainder. *)
      let r = Array.make n 0 in
      for i = 0 to n - 1 do
        r.(i) <- ((u.(i) lsr s) lor (if i + 1 <= n && s > 0 then (u.(i + 1) lsl (bits - s)) land mask else 0))
      done;
      (trim q, trim r)
    end
 
  let to_string a =
    if Array.length a = 0 then "0"
    else
      let rec go a acc = if Array.length a = 0 then acc else let q, r = divmod_small a 10 in
        go q (string_of_int (if Array.length r = 0 then 0 else r.(0)) :: acc) in
      String.concat "" (go a [])
end

Two divisions in base 2^28, both checked against Python's integers, and 100,000 random divisions checked by multiplying back.

10000000000000000000000000000000000000000 / 123456789012345678901
  q = 81000000729000006634
  r = 6661773269766170766
  q*b + r = a: true,  r < b: true
170141183460469231731687303715884105728 / 18446744073709551617
  q = 9223372036854775807
  r = 9223372036854775809
  q*b + r = a: true,  r < b: true
random divisions with q*b + r = a and r < b: 100000 of 100000

The second example, $2^{127} / (2^{64} + 1)$ , is a case in which the estimate survives the refinement and the add-back step runs. Brinch Hansen described how easy it is to get the correction steps subtly wrong, and how rarely the wrong cases occur in testing: an add-back branch that is never exercised can hide a bug for years^[3].

How often the corrections run

In a small base the rare cases become common enough to count. The following measures the first estimate, before refinement, and the add-back rate, over random operands in bases $2^{4}$ , $2^{8}$ and $2^{16}$ :

Statistics for 200,000 random divisions per base.

(* How far the first estimate of each quotient digit is from the true
   digit, and how often the add-back step runs, for small bases where the
   rare cases are frequent enough to count. *)
module type BITS = sig val bits : int end
 
let run (module M : BITS) =
  let module N = D1_lib.Make (M) in
  Random.init 7;
  let rand limbs = N.trim (Array.init limbs (fun _ -> Random.bits () land N.mask)) in
  for _ = 1 to 200_000 do
    let b = rand (2 + Random.int 4) in
    let a = rand (Array.length b + 1 + Random.int 6) in
    if Array.length b >= 2 then ignore (N.divmod a b)
  done;
  let s = N.stats in
  let pct x = 100. *. float x /. float s.digits in
  Printf.printf "base 2^%-2d %8d digits   estimate exact %6.2f%%  1 too large %5.2f%%  2 too large %5.3f%%   add-back %.4f%% (2/base = %.4f%%)\n"
    M.bits s.digits (pct s.qhat_over.(0)) (pct s.qhat_over.(1)) (pct s.qhat_over.(2)) (pct s.add_backs) (200. /. float N.base)

Running it.

base 2^4    871256 digits   estimate exact  75.11%  1 too large 24.48%  2 too large 0.409%   add-back 1.1786% (2/base = 12.5000%)
base 2^8    899177 digits   estimate exact  74.33%  1 too large 24.96%  2 too large 0.714%   add-back 0.0790% (2/base = 0.7812%)
base 2^16   900860 digits   estimate exact  74.39%  1 too large 24.89%  2 too large 0.727%   add-back 0.0003% (2/base = 0.0031%)

The first estimate is exact about three times in four, one too large about a quarter of the time, and two too large in under one per cent of cases, never more, as the bound requires. Knuth estimates the probability that the add-back step is needed as about $2/ β$ ^[1]; with the refinement in place, the measured rate here is about ten times lower, close to $0.2/ β$ in every base. In base $2^{64}$ it is negligible, which is exactly why that branch needs a test constructed to reach it.

§ 05

Complexity and faster methods

Each of the $m + 1$ quotient digits costs one multiply-and-subtract of $n + 1$ limbs, so the algorithm uses $O (mn)$ limb operations, the same as schoolbook multiplication of the quotient by the divisor. For large operands, division is reduced to multiplication. Newton's iteration computes a reciprocal $1/ v$ to the needed precision in a constant number of multiplications, so division costs $O (M (n))$ where $M (n)$ is the cost of multiplying $n$ -limb numbers^[7]. Burnikel and Ziegler's recursive division splits the operands in halves and uses Karatsuba or faster multiplication in the recursive steps, with Algorithm D at the base^[4]. For repeated reduction modulo the same number, Barrett reduction precomputes a scaled reciprocal^[5] and Montgomery multiplication avoids division altogether by working in a transformed representation^[6].

GMP uses schoolbook division, essentially Algorithm D, for small operands and switches to divide-and-conquer and Newton-based methods above thresholds tuned for each processor^[8]. CPython's int implements Algorithm D directly, and Java's BigInteger uses it below a size threshold and Burnikel and Ziegler's method above it.

§ 06

Uses

Every operation in a computer algebra system that reduces a number eventually divides: rational numbers are kept in lowest terms by greatest common divisors, and the Euclidean algorithm is a chain of divisions. Conversion of a large integer to decimal is repeated division by a power of ten. Modular arithmetic for cryptography, and the modular algorithms of computer algebra, need a remainder operation that is exact on every input, including the rare ones that trigger the add-back step. Multiplication above schoolbook size uses Karatsuba or FFT-based methods, and fast division inherits their speed through the reduction to multiplication.

§ 07

History

Long division with estimated digits was part of the earliest multiple-precision arithmetic programs. Pope and Stein studied the accuracy of the digit estimate and the effect of scaling the divisor in 1960^[2]. Knuth gave the algorithm, with the bound on the estimate and the refinement test, as Algorithm D in Seminumerical Algorithms in 1969^[1]. Brinch Hansen's 1994 paper retraced the derivation and documented the pitfalls in implementing it^[3], and Burnikel and Ziegler's 1998 recursive division made it the base case of a subquadratic method^[4].