Here’s bit tricks for 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:
1
2
3
4
5
6
7
8
9
10
11
12
|
uint32_t x = ...
// common method
if ((x%5) == 0)
{
// do something
}
// cheaper method, completely equivalent
if (x * 3435973837U <= 858993459U)
{
\\ do something
}
|
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.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
|
| d |16-bit p|16-bit q| 32-bit p | 32-bit q | 64-bit p | 64-bit q |
|----:|--------|--------|----------:|----------:|--------------------:|--------------------:|
| 3| 43691 | 21845 | 2863311531| 1431655765| 12297829382473034411| 6148914691236517205|
| 5| 52429 | 13107 | 3435973837| 858993459| 14757395258967641293| 3689348814741910323|
| 7| 28087 | 9362 | 3067833783| 613566756| 7905747460161236407| 2635249153387078802|
| 9| 36409 | 7281 | 954437177| 477218588| 10248191152060862009| 2049638230412172401|
| 11| 35747 | 5957 | 3123612579| 390451572| 3353953467947191203| 1676976733973595601|
| 13| 20165 | 5041 | 3303820997| 330382099| 5675921253449092805| 1418980313362273201|
| 15| 61167 | 4369 | 4008636143| 286331153| 17216961135462248175| 1229782938247303441|
| 17| 61681 | 3855 | 4042322161| 252645135| 17361641481138401521| 1085102592571150095|
| 19| 51739 | 3449 | 678152731| 226050910| 9708812670373448219| 970881267037344821|
| 21| 53053 | 3120 | 1022611261| 204522252| 14933078535860113213| 878416384462359600|
| 23| 14247 | 2849 | 3921491879| 186737708| 15238614669586151335| 802032351030850070|
| 25| 23593 | 2621 | 3264175145| 171798691| 10330176681277348905| 737869762948382064|
| 27| 55827 | 2427 | 1749801491| 159072862| 9564978408590137875| 683212743470724133|
| 29| 49717 | 2259 | 1332920885| 148102320| 3816567739388183093| 636094623231363848|
| 31| 31711 | 2114 | 3186588639| 138547332| 17256631552825064415| 595056260442243600|
| 33| 33761 | 1985 | 1041204193| 130150524| 1117984489315730401| 558992244657865200|
| 35| 44939 | 1872 | 2331553675| 122713351| 12649195936257978251| 527049830677415760|
| 37| 7085 | 1771 | 2437684141| 116080197| 1495681951922396077| 498560650640798692|
| 39| 28567 | 1680 | 2532929431| 110127366| 8040888442386214807| 472993437787424400|
| 41| 39961 | 1598 | 3247414297| 104755299| 10348173504763894809| 449920587163647600|
| 43| 48771 | 1524 | 799063683| 99882960| 9437869060967677571| 428994048225803525|
| 45| 20389 | 1456 | 2767867813| 95443717| 5738987045154082725| 409927646082434480|
| 47| 18127 | 1394 | 1736263375| 91382282| 5887258746928580303| 392483916461905353|
| 49| 22737 | 1337 | 438261969| 87652393| 9035139954469984465| 376464164769582686|
| 51| 64251 | 1285 | 4210752251| 84215045| 18085043209519168251| 361700864190383365|
| 53| 21021 | 1236 | 2350076445| 81037118| 2436362424829563421| 348051774975651917|
| 55| 46471 | 1191 | 1483715975| 78090314| 8049488323073258887| 335395346794719120|
| 57| 60937 | 1149 | 3089362441| 75350303| 9385185581360999945| 323627089012448273|
| 59| 55539 | 1110 | 2693454067| 72796055| 14694863923124558067| 312656679215416129|
| 61| 38677 | 1074 | 3238827797| 70409299| 5745707170499696405| 302405640552615600|
| 63| 61375 | 1040 | 3204181951| 68174084| 17275522227759738815| 292805461487453200|
| 65| 4033 | 1008 | 3237744577| 66076419| 1135184250689818561| 283796062672454640|
| 67| 19563 | 978 | 128207979| 64103989| 17345445920055250027| 275324538413575397|
| 69| 4749 | 949 | 2738819725| 62245902| 17377367605668418189| 267344117010283356|
| 71| 43383 | 923 | 3811027319| 60492497| 1818693077689674103| 259813296812810586|
| 73| 61945 | 897 | 3353604601| 58835168| 9097024474706080249| 252695124297391118|
| 75| 51555 | 873 | 2519714147| 57266230| 3443392227092449635| 245956587649460688|
| 77| 14469 | 851 | 1059797125| 55778796| 5749634516480899205| 239568104853370800|
| 79| 5807 | 829 | 1631000239| 54366674| 11208148297950107311| 233503089540627235|
| 81| 18609 | 809 | 2014922929| 53024287| 3188326136196712625| 227737581156908044|
| 83| 17371 | 789 | 724452315| 51746593| 11779246215742243803| 222249928598910260|
| 85| 64765 | 771 | 4244438269| 50529027| 18229723555195321597| 217020518514230019|
| 87| 60263 | 753 | 1875962727| 49367440| 7421103937699244903| 212031541077121282|
| 89| 18409 | 736 | 4198451177| 48258059| 17617676924329347049| 207266787345051141|
| 91| 12243 | 720 | 3539808211| 47197442| 3446095046736949203| 202711473337467600|
| 93| 54261 | 704 | 1062196213| 46182444| 5752210517608354805| 198352086814081200|
| 95| 23455 | 689 | 3571604383| 45210182| 5631111348816599967| 194176253407468964|
| 97| 41889 | 675 | 1594008481| 44278013| 11790702397628785569| 190172619316593315|
| 99| 33099 | 661 | 3210379595| 43383508| 12670490878911611211| 186330748219288400|
| 101| 45421 | 648 | 2083697005| 42524428| 4200743699953660269| 182641030432767837|
|
This table is generated by the Mathematica code
1
2
|
divPair[d_, n_] := {Mod[ExtendedGCD[d, 2^n][[2, 1]], 2^n],Floor[(2^n)/d]};
Table[Flatten[{n, divPair[n, 16], divPair[n, 32], divPair[n, 64]}], {n, 3, 101, 2}]
|
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
Comments
Comment is disabled to avoid unwanted discussions from 'localhost:1313' on your Disqus account...