Wednesday, 6 October 2010

BSDNT - v0.14 divrem_classical

It's time to implement our schoolboy division routine. I prefer the name
divrem_classical, in line with the multiplication routines.

This function will take a dividend a, divisor d, carry-in ci and returns
a quotient q and remainder r. We'll also need to pass in a precomputed
inverse to speed up the dword_t by word_t divisions that need to occur.

We are going to implement the first algorithm spoken of in the previous
post, namely the one which uses a single "normalised" word of the divisor.

We need to be careful not to just shift the top word of the divisor so
that it is normalised, but if it has more than one word, shift any high
bits of the second word as well. We want the top WORD_BITS bits of the
divisor, starting with the first nonzero bit.

We can use our macro for doing a dword_t by word_t division to get each
new word of the quotient. We start with the carry-in and the most
significant word of a. The macro will automatically shift these by the
same amount as we shifted the leading words of the divisor.

As per the divrem1 function, we require that the carry-in be such that
the quotient won't overflow. In other words, we assume that if the
divisor d is m words, then the top m words of the dividend including the
carry-in, are reduced mod d.

First we need to check that we are not in the special case where truncating
the dividend and divisor to two and one words respectively would cause
an overflow of the quotient word to be computed. This only happens
when the top word of the dividend equals the top word of the divisor, as
explained in the previous post.

If the truncation would cause an overflow in the quotient, we collect a
quotient word of ~0, as discussed in the previous post. If not, we compute
the quotient using our macro.

After this point, the remainder is computed. We allow the function to
destroy the input a for this purpose. We leave it up to the caller to make
a copy of a and pass it to the function, if this is not desired
behaviour.

We do our adjustment so that the quotient is correct. We need a while loop
for this, as mentioned in the previous article.

Finally we write out the quotient word and read in the next word of the
dividend that remains.

In the test code we use our mul_classical and muladd_classical functions to
check that divrem_classical is indeed the inverse of these functions, with
zero remainder and nonzero remainder respectively.

The code for today's post is here: v0.14

Previous article: divrem discussion