Programming Praxis

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/(v₁−1)⌋, which is one of the algorithmic tricks used by Knuth to make the code work. (Here, v₁ 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/d₀⌋, 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