## Big Numbers: Division

### June 7, 2011

We’ll show only the main function that performs division of two positive integers; you can read the rest of the function on the next page. Here it is:

```(define (div ns ds)   (if (lt? ns ds) (values '(0) ns)     (let* ((n (length ds)) (m (- (length ns) n)))       (if (= n 1) (div1 ns (car ds))         (let* ((d (quotient big-base (+ (car (reverse ds)) 1)))                (rev-ns (reverse (mul1 ns d)))                (rev-ns (if (= (length ns) (length rev-ns)) (cons 0 rev-ns) rev-ns))                (ds (mul1 ds d)) (d0 (car (reverse ds))) (d1 (cadr (reverse ds))))           (let loop ((j m) (rev-ns rev-ns) (qs '()))             (if (negative? j)                 (call-with-values                   (lambda () (div1 (reverse rev-ns) d))                   (lambda (q r) (values (reverse (reduce (reverse qs))) q)))                 (let* ((x (+ (* (car rev-ns) big-base) (cadr rev-ns)))                        (q (nextq x d0 d1 rev-ns))                        (ds*q (mul1 ds q)))                   (if (lt? (reverse (take (+ n 1) rev-ns)) ds*q)                       (let* ((q (- q 1)) (ds*q (sub ds*q ds)))                         (loop (- j 1) (nextn j n rev-ns ds*q) (cons q qs)))                       (loop (- j 1) (nextn j n rev-ns ds*q) (cons q qs)))))))))))```

First we pick off two special cases: where n < d, the quotient is 0 and the remainder is n, and where d has only a single digit, we can use the simpler function `div1` (we need `div1` later, anyway).

The second `let*` block performs normalization, multiplying n and d by ⌊b/(v1−1)⌋, which is one of the algorithmic tricks used by Knuth to make the code work. (Here, v1 is Knuth’s notation for the most significant digit in the divisor. Confusingly, we use d for the divisor and Knuth uses d for the coefficient of normalization. Sorry!)

The main loop tracks three quantities: `j` is a count of the digits in the quotient, `rev-ns` is the remaining digits in the partial dividend, in reverse order, and `qs` is the accumulating digits of the quotient. The loop ends when `j` becomes negative, returning the quotient and the denormalized partial dividend as the remainder.

The main body of the loop computes the next digit of the quotient at each step. The `let*` computes `x`, which is the first two digits of the partial dividend, expressed as a number (it must always be less than the square of the radix base), `q`, which is the proposed next digit of the quotient, and `ds*q`, which is the divisor times the proposed next digit of the quotient. In Knuth’s algorithm, `q` (he calls it qhat) is special; it is either the true next digit of the quotient, or at most too large by one. Here is the computation of the next `q`:

```(define (nextq x d0 d1 ns)   (let loop ((q (quotient x d0)))     (if (< (* q d1) (+ (* (- x (* q d0)) big-base) (caddr ns))) q       (loop (- q 1)))))```

Here `x` is the first two digits of the partial dividend expressed as a number, `d0` and `d1` are the two most significant digits of the divisor, and `ns` is the remaining partial dividend. `Q` is first estimated as ⌊x/d0⌋, then repeatedly reduced by 1 as long as the first three digits of the partial dividend are greater than `q` times the first two digits of the divisor. Knuth proved that this calculation always makes `q` greater than or equal to the true next digit of the quotient (Theorem 4.3.1 A) but never leaves `q` too large by more than one (Theorem 4.3.1 B). The clever calculation also assures that no intermediate quantities are ever larger than the square of the radix base.

The next calculation in the main loop compares `q` times the full divisor (not just the first two digits) to the full partial dividend (not just the first three digits), and reduces `q` by 1 if necessary. Then the code loops. The tricky part here is the calculation of the new value of the partial dividend, which may have to add a leading 0 to the partial dividend to preserve the sense of the first three digits:

```(define (nextn j n rev-ns ds*q)   (let ((zs (append (reverse (sub (reverse (take (+ n 1) rev-ns)) ds*q))                     (drop (+ n 1) rev-ns))))     (if (< (length zs) (+ n j)) (cons 0 zs) zs)))```

And that’s it. The full details are on the next page, where you can see the basic addition, subtraction and multiplication code, as well as `mul1` and `div1`, which take a big number and a single big-digit. That code also shows the error generated on division by zero and the calculation of the signs of the quotient and remainder when either the dividend or divisor are negative.

In addition to `big-divide`, we also provide `big-quotient` and `big-remainder`, by analogy to Scheme:

```(define (big-quotient big1 big2)   (call-with-values     (lambda () (big-divide big1 big2))     (lambda q r) q))```

```(define (big-remainder big1 big2)   (call-with-values     (lambda () (big-divide big1 big2))     (lambda q r) r))```

You can run the program at http://programmingpraxis.codepad.org/MJl0kPSU. The radix was set to 10 and debugging code was turned on so you can see all the details of how the function works.

Pages: 1 2 3

### 2 Responses to “Big Numbers: Division”

1. […] today’s Programming Praxis exercise, our goal is to add division to our Big Number library. Let’s get […]

```import Control.Arrow

instance Integral BigNum where
quotRem n@(B l1 ds1) d@(B l2 ds2)
| d == 0         = error "Division by zero"
| n < 0 || d < 0 = let (q',r') = quotRem (abs n) (abs d)
in (signum n * signum d * q', signum n * r')
| n < d          = (0, n)
| otherwise      = first (+ q) \$ quotRem (n - d*q) d where
(n1,n2,d1) = (last ds1, last \$ tail ds1, last ds2)
q = if n1 < d1 then shift (l1 - l2 - 1) \$ div n2 d1 +
fromIntegral (div base (fromIntegral d1))
else shift (l1 - l2) \$ div n1 d1
shift s i = B (s + 1) \$ (replicate s 0) ++ [i]
```