Saturday, 18 September 2010

BSDNT - v0.9 divrem1_hensel

The next division function we'll introduce is so-called Hensel division, or
right-to-left division. This starts at the least significant word and applies
the usual division algorithm back towards the most significant word.

Of course, now any "remainder" will be at the most significant word end and
have a completely different meaning.

Hensel division is actually just division mod B^m. Thus if q is the Hensel
quotient, you have a = q * d mod B^m. In fact, if r is the Hensel remainder,
you have a = q*d + r*B^m.

If a fits in m words and the division is exact (i.e. a = q*d), the Hensel
division and ordinary (euclidean) division return the same quotient. If the
division is not exact, this is no longer the case.

However, it may be that chaining Hensel divisions together does produce an
exact division, even when division over the bottom part of the chain
wouldn't produce an exact division. Thus, the ability to chain Hensel
divisions, and thus the ability to return the Hensel remainder, is

Well, now we come to the first issue if we want to implement this. What is
the C operator for Hensel division? Actually, there doesn't seem to be one.
We have to implement it ourselves.

If a0 is a word from our dividend, and d is our single word divisor, then
we require q0 such that d * q0 = a0 mod B. In other words, we want q0 and
r0 such that d * q0 = a0 + r0*B.

Suppose we can find a value dinv such that dinv * d = 1 mod B. Then
we have that dinv * d * q0 = dinv * a0 mod B, i.e. q0 = dinv * a0 mod B,
which is a single word-by-word mullow.

Once we have q0 we can compute d * q0, the high word of which is the Hensel
"remainder" r0. We need to subtract this from the next limb of our divisor
before continuing.

We'll make the convention that carry-ins and carry-outs from our Hensel
division are positive rather than negative. We'll just remember to subtract
them. This way, if q is our Hensel quotient and r our Hensel remainder,
then d * q = a + r * B^n.

So how do we compute dinv? Firstly, we better require that d is odd, else
dinv * d = 1 mod B is an impossibility. With this restriction, we have
gcd(d, B) = 1 and therefore the extended euclidean algorithm lets us find
s, t such that s*d + t*B = 1. Then modulo B we have s*d = 1 mod B, i.e.
dinv = s is the precomputed inverse that we are after.

Another way to solve dinv * d = 1 mod B is to use Hensel lifting. This
works by first solving the equation mod 2, then use that solution to
solve it mod 4 and so on. This can all be done with nothing more than

Firstly, we know a solution mod 2, namely 1 (as d is odd), i.e. if v1 = 1
then we have v1 * d = 1 mod 2. Now suppose we have a solution vk mod 2^k,
i.e. vk * d = 1 mod 2^k. We'd like to find a value wk such that
(vk + 2^k * wk) * d = 1 mod 2^2k. This is always possible -- I omit the
easy proof -- so we only need to solve this equation for wk, i.e.
(d * vk - 1) = -2^k * wk * d. We can get rid of the d on the right hand
side by multiplying by vk, i.e. 2^k * w = ((1 - d * vk) * vk) (mod 2^2k).

In other words, two multiplications will suffice to compute
2^k * wk (mod 2^2k) from vk. Then v{k+1} = vk + 2^k * wk is the inverse
of d mod 2^2k.

As a further optimisation of all this, one can just work mod 2^64 all
the way through the computation. In other words, start with v1 = 1
and compute v{k+1} = vk + ((1 - d * vk) * vk) (mod 2^64), six times.

To get from a solution mod 2 to a solution mod 2^64 clearly only takes
six Hensel lifting steps. Of course, with a lookup table mod 2^8 to start
from, instead of starting from a solution mod 2, one can do it in three
steps, requiring just six word-by-word multiplications (mullow only).

However, today I am feeling very lazy and I leave it as an exercise for
someone to implement the Hensel lifting lookup table method.

We add a simple type to hold our precomputed Hensel inverse, and a function
for computing it.

The actual nn_divrem_hensel_1_preinv function is again a straightforward
loop. We ensure our remainder is positive at each step and ensure that
we propagate any borrows from subtractions at each step.

The test code for the Hensel division is very similar to that for the
ordinary euclidean division except that d must be odd and the chaining
test will proceed right-to-left.

OK, that is enough brain strain for one day!

Here v0.9 is the github repository for today's post.

Previous post: v0.8 - divrem1_preinv

1 comment: