Efficient divisibility testing

Efficient divisibility testing

The example

Here is a low computational cost method to test if an integer $x$ is divisible by another integer $d$, useful in C-like languages (assembly included). It replaces a possibly costly modulus % operation with a usually lower cost multiplication *.

Here is an example illustrating the method for divisibility by 5:

What trickery is this? How does it work? Where do the constants 3435973837U and 858993459U come from?

Note the magic numbers depend on the type of $x$. The above example requires $x$ is an unsigned 32 bit integer, and the language has multiplication overflow computed modulus 32. Many common languages (C,C++,C#,Java, ..) do this.

The values

Before giving the details, here is a table of the magic constants for divisibility testing of odd integers from 3 to 101 inclusive. The bitsize is the size of the unsigned integer $x$, and the test is then x*p <= q.

This table is generated by the Mathematica code

The method

Here is how it works, and how you can compute your own values.

Suppose $x$ is an $n$-bit unsigned integer, and write $N=2^n$. Then for any odd integer divisor $d < N$ there is always an integer $p$ depending on $d$ (called a multiplicative inverse mod $N$) so that $pd$ is one more than a multiple of $N$, i.e., $pd=1 + N k$ for some positive integer $k$. This $p$ can be found via the Extended Euclidean algorithm for computing a GCD [1].

The value $q$ to compare against is computed as $q=\lfloor N/d \rfloor$, where $\lfloor y \rfloor$ is the floor function (largest integer less than or equal to $y$).

Ok, so those are the values to use. How does it work? Here is the rough idea, then a proof.

The idea

In the ring of integers mod $N$, multiplying by $p$ acts like dividing by $d$, since for any integer $x$, $(xp)d\equiv x(pd)\equiv x (1)\equiv x$ mod $N$. So think of $xp$ as $x/d$ mod $N$.

The comparison value is the value $N/d$, taken again in the ring of integers mod $N$. The multiplication by $p$ "shuffles" the integers in the ring, so those divisible by $d$ are ordered first, then those with "residue" 1, then those with residue 2, etc. It’s weird, but it works.

The proof

Here is a proof of why it works. Write $\left[t\right]$ to treat the integer $t$ mod $N$ as an integer again for comparison purposes. We need to show $d$ divides $x$ if and only if $[x\cdot p]\leq q$, with $p$ and $q$ as above.

Let $b$ be the remainder of $xp$ divided by $N$, i.e., there is an integer $a$ and an integer integer $0\leq b<N$, with $xp=Na+b$. Thus $[xp] = b$.

Assume $d$ divides $x$ (i.e., $x= t d$ for some integer $t$). Then $[xp]=[tdp]=[t]=t$, and $x<N$ implies $t<N/d$ implies $t\leq q = \lfloor N/d\rfloor$. Thus $[xp]\leq q$ as required.

Conversely, suppose $[xp]\leq q$. Then $xp=Na+b$ implies (expanding two ways) $xpd = x(1+Nk)=x+Nkx$ and $xpd=(Na+b)d =bd + Nad$. Thus $x+Nkx = bd + Nad$. But $b=[xp]\leq q$ implies $bd\leq qd = \lfloor N/d/rfloor d < N$. The strict inequality follows from $N$ and $d$ being relatively prime. Thus, since $0\leq bd < N$, taking both sides of the double expansion mod $N$, we get $x=bd$.

Q.E.D.

The end

This can be extended to divisibility testing against any $d$ (odd or even) by first testing $x$ against the largest power of 2 dividing $d$, then shifting that out if it passed and testing the odd part.

[1] https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

  • Hans Bezemer

    You’re a hard man to find 😉 I wanted to thank you for your WELL512 implementation, which I translated to Forth, in particular: my Forth implementation “4tH”. It will be released this year as v3.62.5. If you’re interested, you can google “4tH compiler” or “4tH compiler sourceforge” if you’re particularly interested to see your code again 😉

  • Chris Lomont

    Cool – glad you liked it. If you like random number generation topics, I just posted another one at http://clomont.com/generating-uniform-random-numbers-on-01/