# Lomont.org

(Also lomonster.com and clomont.com)

Published

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.

Categories:

Tags:

## Comments

Comment is disabled to avoid unwanted discussions from 'localhost:1313' on your Disqus account...