Thursday, 13 September 2012

Divisibility and modular multiplication

An other post from CodeProject (with small changes). I promise I'll post something new next time.

Typically, all computer math is done modulo 232 (or 2something, which the rest of my post trivially generalizes to). This leads to sometimes surprising effects, such as that the sum of `a` and `b`, while both positive, can end up being lower than `min(a, b)`, which is rather well-known and known as "overflow" (often treated like something bad, which it can be). Less well known is that it also means that some numbers have a multiplicative inverse, ie a number `x-1` such that `x-1x=1`. As mentioned in the wikipedia article, the numbers which have multiplicative inverses are precisely those coprime to the modulo. The modulo is a power of two, which means that a number has a multiplicative inverse iff it is odd.

And it turns out to be actually useful, too. One application is, as you might guess from the title, divisibility testing.

Multiplying a number by an odd number is reversible - you can take the multiplicative inverse of the odd number and then multiply by it to get the original number back. Put differently, the function `f(x) = x * k (for k odd)` is a bijection.

Modular multiplication is associative, so a multiple of `k`, say `n * k` multiplied by `inv(k)` (the multiplicative inverse of `k`), is `n`, because

```(n * k) * inv(k) =
// use associativity
n * (k * inv(k)) =
// use definition of multiplicative inverse
n * 1 =
// multiplicative identity
n```
That means that in `x * inv(k)`, the multiples of `k` "use up" the results from 0 through (232-1)/`k`, they can't be used twice because it's a bijection, leaving just the numbers bigger than (232-1)/`k` for non-multiples-of-k. Which suggest a very simple divisibility test:
```static bool IsDivisibleByOdd(uint x, uint divisor)
{
if ((divisor & 1) == 0)
throw new ArgumentException("divisor must be odd");
uint d_inv = inv(divisor);
uint biggest = uint.MaxValue / divisor;  // problem right here
return (x * d_inv) <= biggest;
}

static uint inv(uint d)
{
// see Hacker's Delight,
// Computing the Multiplicative Inverse by Newton's Method
// use extra iteration when extending to 64 bits
uint x = (d * d) + d - 1;
uint t = d * x;
x *= 2 - t;
t = d * x;
x *= 2 - t;
t = d * x;
x *= 2 - t;
return x;
}```

This may seem at first glance not to help at all, but for a constant divisor all the scary operations can be precomputed. It even has some use for unknown divisors, as long as the inverse and upper limit can be reused often enough.

I will extend this method to even numbers after the next post, which will cover some concept (bits again, I haven't forgotten about them) that are necessary for it.