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.

New enough versions of GCC and Clang can perform this optimization when the divisor is a constant, but not yet in cases where the same (but unknown) divisor is re-used often.

This method can be extended to even divisors, with some complications.