Unnormalizing in Knuth's Algorithm D

Question

I'm trying to implement Algorithm D from Knuth's "The Art of Computer Programming, Vol 2" in Rust although I'm having trouble understating how to implement the very last step of unnormalizing. My natural numbers are a class where each number is a vector of u64, in base u64::MAX. Addition, subtraction, and multiplication have been implemented.

Knuth's Algorithm D is a euclidean division algorithm which takes two natural numbers x and y and returns (q,r) where q = x / y (integer division) and r = x % y, the remainder. The algorithm depends on an approximation method which only works if the first digit of y is greater than b/2, where b is the base you're representing the numbers in. Since not all numbers are of this form, it uses a "normalizing trick", for example (if we were in base 10) instead of doing 200 / 23, we calculate a normalizer d and do (200 * d) / (23 * d) so that 23 * d has a first digit greater than b/2.

So when we use the approximation method, we end up with the desired q but the remainder is multiplied by a factor of d. So the last step is to divide r by d so that we can get the q and r we want. My problem is, I'm a bit confused at how we're suppose to do this last step as it requires division and the method it's part of is trying to implement division.

(Maybe helpful?): The way that d is calculated is just by taking the integer floor of b-1 divided by the first digit of y. However, Knuth suggests that it's possible to make d a power of 2, as long as d * the first digit of y is greater than b / 2. I think he makes this suggestion so that instead of dividing, we can just do a binary shift for this last step. Although I don't think I can do that given that my numbers are represented as vectors of u64 values, instead of binary.

Any suggestions?

Answer 1

You are right, you should pick d as a power of two. You need just a shift: shift = count_leading_zeroes(divisor[divisor.size() - 1]). Then, instead of multiplying and dividing, you can shift your big integers. Here is an example code for shifting big integers:

// For base = 2^64 and 0 <= shift <= 64
void shift_right_inplace(bint_t& a, int shift) {
    for (int i = 0; i < a.data.size(); i++) {
        const uint64_t self = a.data[i] >> shift;
        const uint64_t prev = i == a.data.size() - 1 ? 0 : a.data[i + 1] << (64 - shift);
        a.data[i] = self | prev;
    }
    if (a.data.back() == 0) {
        a.data.pop_back();
    }
}

For example, let's take base as 256 for simplicity and look at two adjacent "digits": abcdefgh and klmnopqr. You want to shift it right by 3.

The first byte has no neighbor from the left, so you just shift it: 000abcde. For the second byte, you shift it (000klmno) and take the rest from the previous byte (abcdefgh >> 3 = 000abcde | fgh, it is effectively prev << (8 - 3)).