Faster Remainder by Direct Computation
\subtitlefontApplications to Compilers and Software Libraries

Integer division often refers to two closely related concepts, the actual division and the modulus. Given an integer numerator $n$ and a non-zero integer divisor $d$, the integer division, written $\operatorname*{\mathrm{div}}$, gives the integer quotient ($n\operatorname*{\mathrm{div}}d=q$). The modulus, written $\bmod$, gives the integer remainder ($n\bmod d=r$). Given an integer numerator $n$ and an integer divisor $d$, the quotient ($q$) and the remainder ($r$) are always integers even when the fraction $n/d$ is not an integer. It always holds that the quotient multiplied by the divisor plus the remainder gives back the numerator: $n=q\operatorname{\ast}d+r$.

Depending on the context, ‘integer division’ might refer solely to the computation of the quotient, but might also refer to the computation of both the integer quotient and the remainder. The integer division instructions on x64 processors compute both the quotient and the remainder.^aaaWe use x64 to refer to the commodity Intel and AMD processors supporting the 64-bit version of the x86 instruction set. It is also known as x86-64, x86_64, AMD64 and Intel 64. In most programming languages, they are distinct operations: the C programming language uses / for division ($\operatorname*{\mathrm{div}}$) and % for modulo ($\bmod$).

Let us work through a simple example to illustrate how we can replace an integer division by a multiplication. Assume we have a pile of 23 items, and we want to know how many piles of 4 items we can divide it into and how many will be left over ($n=23$, $d=4$; find $q$ and $r$). Working in base 10, we can calculate the quotient $23\operatorname*{\mathrm{div}}4=5$ and the remainder $23\bmod 4=3$, which means that there will be 5 complete piles of 4 items with 3 items left over ($23=5\operatorname{\ast}4+3$).

If for some reason, we do not have a runtime integer division operator (or if it is too expensive for our purpose), we can instead precompute the multiplicative inverse of 4 once ($c=1/d=1/4=0.25$) and then calculate the same result using a multiplication ($23\operatorname{\ast}0.25=5.75$). The quotient is the integer portion of the product to the left of the decimal point ($q=\left\lfloor 5.75\right\rfloor=5$), and the remainder can be obtained by multiplying the fractional portion $f=0.75$ of the product by the divisor $d$: $r=f\operatorname{\ast}d=0.75\operatorname{\ast}4=3$.

The binary registers in our computers do not have a built-in concept of a fractional portion, but we can adopt a fixed-point convention. Assume we have chosen a convention where $1/d$ has 5 bits of whole integer value and 3 bits of ‘fraction’. The numerator 23 and divisor 4 would still be represented as standard 8-bit binary values (00010111 and 00000100, respectively), but $1/d$ would be 00000.010. From the processor’s viewpoint, the rules for arithmetic are still the same as if we did not have a binary point—it is only our interpretation of the units that has changed. Thus we can use the standard (fast) instructions for multiplication ($00010111\operatorname{\ast}00000010=00101110$) and then mentally put the ‘binary point’ in the correct position, which in this case is 3 from the right: 00101.110. The quotient $q$ is the integer portion (leftmost 5 bits) of this result: 00101 in binary ($q=5$ in decimal). In effect, we can compute the quotient $q$ with a multiplication (to get 00101.110) followed by a right shift (by three bits, to get 000101). To find the remainder, we can multiply the fractional portion (rightmost 3 bits) of the result by the divisor: $00000.110\operatorname{\ast}00000100=00011.000$ ($r=3$ in decimal). To quickly check whether a number is divisible by 4 ($n\mod 4=0$?) without computing the remainder it suffices to check whether the fractional portion of the product is zero.

But what if instead of dividing by 4, we wanted to divide by 6? While the reciprocal of 4 can be represented exactly with two digits of fixed-point fraction in both binary and decimal, $1/6$ cannot be exactly represented in either. As a decimal fraction, $1/6$ is equal to the repeating fraction 0.1666…(with a repeating 6), and in binary it is 0.0010101…(with a repeating 01). Can the same technique work if we have a sufficiently close approximation to the reciprocal for any divisor, using enough fractional bits? Yes!

For example, consider a convention where the approximate reciprocal $c$ has 8 bits, all of which are fractional. We can use the value 0.00101011 as our approximate reciprocal of 6. To divide $23$ by 6, we can multiply the numerator (10111 in binary) by the approximate reciprocal: $n\operatorname{\ast}c=00010111\operatorname{\ast}0.00101011=11.11011101$. As before, the decimal point is merely a convention, the computer need only multiply fixed-bit integers. From the product, the quotient of the division is 11 in binary ($q=3$ in decimal); and indeed $23\operatorname*{\mathrm{div}}6=3$. To get the remainder, we multiply the fractional portion of the product by the divisor ($f\operatorname{\ast}d=0.11011101\operatorname{\ast}00000110=101.00101110$), and then right shift by 8 bits, to get 101 in binary ($r=5$ in decimal). See Table 1 for other examples.

While the use of the approximate reciprocal $c$ prevents us from confirming divisibility by 6 by checking whether the fractional portion is exactly zero, we can still quickly determine whether a number is divisible by 6 ($n\bmod 6=0$?) by checking whether the fractional portion is less than the approximate reciprocal ($f<c$?). Indeed, if $n=q\operatorname{\ast}d+r$ then the fractional portion of the product of $n$ with the approximate reciprocal should be close to $r/d$: it makes intuitive sense that comparing $r/d$ with $c\approx 1/d$ determines whether the remainder $r$ is zero. For example, consider $n=42$ (101010 in binary). We have that our numerator times the approximate reciprocal of 6 is $101010\operatorname{\ast}0.00101011=111.00001110$. We see that the quotient is 111 in binary ($q=7$ in decimal), while the fractional portion $f$ is smaller than the approximate reciprocal $c$ ($0.00001110<0.00101011$), indicating that 42 is a multiple of 6.

In our example with 6 as the divisor, we used 8 fractional bits. The more fractional bits we use, the larger the numerator we can handle. An insufficiency of fractional bits can lead to incorrect results when $n$ grows. For instance, with $n=131$ ($10000011$ in binary) and only 8 fractional bits, $n$ times the approximate reciprocal of 6 is $10000011\operatorname{\ast}0.00101011=10110.00000001$, or 22 in decimal. Yet in actuality, $131\operatorname*{\mathrm{div}}6=21$; using an 8 bit approximation of the reciprocal was inadequate.

How close does the approximation need to be?—that is, what is the minimum number of fractional bits needed for the approximate reciprocal $c$ such that the remainder is exactly correct for all numerators? We derive the answer in § 3.

The scenario we describe with an expensive division applies to current processors. Indeed, integer division instructions on recent x64 processors have a latency of 26 cycles for 32-bit registers and at least 35 cycles for 64-bit registers ¹. We find similar latencies in the popular ARM processors ². Thus, most optimizing compilers replace integer divisions by constants $d$ that are known at compile time with the equivalent of a multiplication by a scaled approximate reciprocal $c$ followed by a shift. To compute the remainder by a constant ($n\bmod d$), an optimizing compiler might first compute the quotient $n\operatorname*{\mathrm{div}}d$ as a multiplication by $c$ followed by a logical shift by $F$ bits $(c\operatorname{\ast}n)\operatorname*{\mathrm{div}}2^{F}$, and then use the fact that the remainder can be derived using a multiplication and a subtraction as $n\bmod d=n-(n\operatorname*{\mathrm{div}}d)\operatorname{\ast}d$.

Current optimizing compilers discard the fractional portion of the multiplication ($(c\operatorname{\ast}n)\bmod 2^{F}$). Yet using the fractional bits to compute the remainder or test the divisibility in software has merit. It can be faster (e.g., by more than 25%) to compute the remainder using the fractional bits compared to the code produced for some processors (e.g., x64 processors) by a state-of-the-art optimizing compiler.

Jacobsohn³ derives an integer division method for unsigned integers by constant divisors. After observing that any integer divisor can be written as an odd integer multiplied by a power of two, he focuses on the division by an odd divisor. He finds that we can divide by an odd integer by multiplying by a fractional inverse, followed by some rounding. He presents an implementation solely with full adders, suitable for hardware implementations. He observes that we can get the remainder from the fractional portion with rounding, but he does not specify the necessary rounding or the number of bits needed.

In the special case where we divide by 10, Vowels⁴ describes the computation of both the quotient and remainder. In contrast with other related work, Vowels presents the computation of the remainder directly from the fractional portion. Multiplications are replaced by additions and shifts. He does not extend the work beyond the division by 10.

Granlund and Montgomery⁵ present the first general-purpose algorithms to divide unsigned and signed integers by constants. Their approach relies on a multiplication followed by a division by a power of two which is implemented as an logical shift ($(\left\lceil 2^{F}/d\right\rceil\operatorname{\ast}n)\operatorname*{\mathrm{div}}2^{F}$). They implemented their approach in the GNU GCC compiler, where it can still be found today (e.g., up to GCC version 7). Given any non-zero 32-bit divisor known at compile time, the optimizing compiler can (and usually will) replace the division by a multiplication followed by a shift. Following Cavagnino and Werbrouck⁶, Warren⁷ finds that Granlund and Montgomery choose a slightly suboptimal number of fractional bits for some divisors. Warren’s slightly better approach is found in LLVM’s Clang compiler. See Algorithm 1.

Following Artzy et al.⁸, Granlund and Montgomery⁵ describe how to check that an integer is a multiple of a constant divisor more cheaply than by the computation of the remainder. Yet, to our knowledge, no compiler uses this optimization. Instead, all compilers that we tested compute the remainder $r$ by a constant using the formula $r=n-q\operatorname{\ast}d$ and then compare against zero. That is, they use a constant to compute the quotient, multiply the quotient by the original divisor, subtract from the original numerator, and only finally check whether the remainder is zero.

In support of this approach, Granlund and Montgomery⁵ state that the remainder, if desired, can be computed by an additional multiplication and subtraction. Warren⁷ covers the computation of the remainder without computing the quotient, but only for divisors that are a power of two $2^{K}$, or for small divisors that are nearly a power of two ($2^{K}+1$, $2^{K}-1$).

In software, to our knowledge, no authors except Jacobsohn³ and Vowels⁴ described using the fractional portion to compute the remainder or test the divisibility, and neither of these considered the general case. In contrast, the computation of the remainder directly, without first computing the quotient, has received some attention in the hardware and circuit literature^{9, 10, 11}. Moreover, many researchers^{12, 13} consider the computation of the remainder of unsigned division by a small divisor to be useful when working with big integers (e.g., in cryptography).

Instead of discarding the least significant bits resulting from the multiplication in the Granlund-Montgomery-Warren algorithm, we can use them to compute the remainder without ever computing the quotient. We formalize this observation by Theorem 3.3, which we believe to be a novel mathematical result.

In general, we expect that it takes at least $N$ fractional bits for the approximate reciprocal $c$ to provide exact computation of the remainder for all non-negative numerators less than $2^{N}$. Let us say we use $F=N+L$ fractional bits for some non-negative integer value $L$ to be determined. We want to pick $L$ so that approximate reciprocal $c=\left\lceil 2^{F}/d\right\rceil=\left\lceil 2^{N+L}/d\right\rceil$ allows exact computation of the remainder as $r=\left(\left(\left(c\operatorname{\ast}n\right)\bmod 2^{F}\right)\operatorname{\ast}d\right)\operatorname*{\mathrm{div}}2^{F}$ where $F=N+L$.

We illustrate our notation using the division of $63$ by $6$ as in the last row of Table 1. We have that $63$ is 111111 in binary and that the approximate reciprocal $c$ of $6$ is $\left\lceil 2^{8}/6\right\rceil=00101011$ in binary. We can compute the quotient as the integer part of the product of the reciprocal by

Taking the $F$-bit fractional portion ($(c\operatorname{\ast}n)\bmod 2^{F}$), and multiplying it by the divisor $6$ (110 in binary), we get the remainder as the integer portion of the result:

The fractional portion $(c\operatorname{\ast}n)\bmod 2^{F}$ given by 10010101 is relatively close to the product of the reciprocal by the remainder ($00101011\operatorname{\ast}11$) given by 10000001: as we shall see, this is not an accident.

Indeed, we begin by showing that the $F=N+L$ least significant bits of the product ($(c\operatorname{\ast}n)\bmod 2^{N+L}$) are approximately equal to the scaled approximate reciprocal $c$ times the remainder we seek ($c\operatorname{\ast}(n\bmod d)$), in a way made precise by Lemma 3.1. Intuitively, this intermediate result is useful because we only need to multiply this product by $d$ and divide by $2^{F}$ to cancel out $c$ (since $c\operatorname{\ast}d\approx 2^{F}$) and get the remainder $n\bmod d$.

Lemma 3.1 tells us that $(c\operatorname{\ast}n)\bmod 2^{N+L}$ is close to $c\operatorname{\ast}(n\bmod d)$ when $c$ is close to $2^{N+L}/d$. Thus it should make intuitive sense that $(c\operatorname{\ast}n)\bmod 2^{N+L}$ multiplied by $d/2^{N+L}$ should give us $n\bmod d$. The following theorem makes the result precise.

Consider the constraint $2^{N+L}\leq c\operatorname{\ast}d\leq 2^{N+L}+2^{L}$ given by Theorem 3.3.

• •

  We have that $c=\left\lceil 2^{N+L}/d\right\rceil$ is the smallest value of $c$ satisfying $2^{N+L}\leq c\operatorname{\ast}d$.

• •

  Furthermore, when $d$ does not divide $2^{N+L}$, we have that $\left\lceil 2^{N+L}/d\right\rceil\operatorname{\ast}d=2^{N+L}+d-(2^{N+L}\bmod d)$ and so $c\operatorname{\ast}d\leq 2^{N+L}+2^{L}$ implies $d\leq(2^{N+L}\bmod d)+2^{L}$. See Lemma 3.5.

On this basis, Algorithm 2 gives the minimal number of fractional bits $F$. It is sufficient to pick $F\geq N+\log_{\lx@text@underscore}2(d)$.

It is not always best to use the smallest number of fractional bits. For example, we can always conveniently pick $F=2\operatorname{\ast}N$ (meaning $L=N$) and $c=\left\lceil 2^{2N}/d\right\rceil$, since $d\leq 2^{F}\bmod d+2^{N}$ clearly holds (given $d\leq 2^{N}$).