Faster Remainder by Direct Computation: Applications to Compilers and Software Libraries

TL;DR

This paper introduces a new algorithm for directly computing remainders using inverse multiplication, leading to faster division and divisibility tests on modern processors, outperforming existing compiler-optimized code.

Contribution

A novel, generally applicable method for directly computing remainders and divisibility tests that surpasses current compiler optimizations in speed.

Findings

Our method is up to 25% faster than compiler-optimized code.

Divisibility tests are up to 50% faster than state-of-the-art methods.

The approach works efficiently across multiple processor architectures.

Abstract

On common processors, integer multiplication is many times faster than integer division. Dividing a numerator n by a divisor d is mathematically equivalent to multiplication by the inverse of the divisor (n / d = n x 1/d). If the divisor is known in advance---or if repeated integer divisions will be performed with the same divisor---it can be beneficial to substitute a less costly multiplication for an expensive division. Currently, the remainder of the division by a constant is computed from the quotient by a multiplication and a subtraction. But if just the remainder is desired and the quotient is unneeded, this may be suboptimal. We present a generally applicable algorithm to compute the remainder more directly. Specifically, we use the fractional portion of the product of the numerator and the inverse of the divisor. On this basis, we also present a new, simpler divisibility algorithm to detect nonzero remainders. We also derive new tight bounds on the precision required when representing the inverse of the divisor. Furthermore, we present simple C implementations that beat the optimized code produced by state-of-art C compilers on recent x64 processors (e.g., Intel Skylake and AMD Ryzen), sometimes by more than 25%. On all tested platforms including 64-bit ARM and POWER8, our divisibility-test functions are faster than state-of-the-art Granlund-Montgomery divisibility-test functions, sometimes by more than 50%.

Figures

Tables

Table 1. Table 1: Division by 6 ( d = 110 𝑑 110 d=110 in binary) using a multiplication by the approximate reciprocal ( c = 0.00101011 𝑐 0.00101011 c=0.00101011 in binary). The numerator n 𝑛 n is an N 𝑁 N -bit value, with N = 6 𝑁 6 N=6 . The approximate reciprocal uses F = 8 𝐹 8 F=8  fractional bits. The integer portion of the product ( N 𝑁 N bits) gives the quotient. Multiplying the fractional portion of the product ( F 𝐹 F bits) by the divisor ( N 𝑁 N bits) and keeping only the integer portion ( N 𝑁 N bits), we get the remainder (two last columns). The integer portion in bold (column 2) is equal to the quotient (column 3) in binary. The integer portion in bold (column 4) is equal to the remainder (column 5) in binary.

$n$	numerator times the approx. reciprocal ($n\ast c$)	quotient	fractional portion $\ast$ divisor ($f\ast d$)	remainder
	$N$ bits $\ast$ $F$ bits $\to$ $N+F$ bits	$N$ bits	$F$ bits $\ast$ $N$ bits $\to$ $N+F$ bits	$N$ bits
0	$000000\ast 0.00101011=\text{000000}.00000000$	0	$0.00000000\ast 000110=\text{000000}.00000000$	0
1	$000001\ast 0.00101011=\text{000000}.00101011$	0	$0.00101011\ast 000110=\text{000001}.00000010$	1
2	$000010\ast 0.00101011=\text{000000}.01010110$	0	$0.01010110\ast 000110=\text{000010}.00000100$	2
3	$000011\ast 0.00101011=\text{000000}.10000001$	0	$0.10000001\ast 000110=\text{000011}.00000110$	3
4	$000100\ast 0.00101011=\text{000000}.10101100$	0	$0.10101100\ast 000110=\text{000100}.00001000$	4
5	$000101\ast 0.00101011=\text{000000}.11010111$	0	$0.11010111\ast 000110=\text{000101}.00001010$	5
6	$000110\ast 0.00101011=\text{000001}.00000010$	1	$0.00000010\ast 000110=\text{000000}.00001100$	0
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
17	$010001\ast 0.00101011=\text{000010}.11011011$	2	$0.11011011\ast 000110=\text{000101}.00100010$	5
18	$010010\ast 0.00101011=\text{000011}.00000110$	3	$0.00000110\ast 000110=\text{000000}.00100100$	0
19	$010011\ast 0.00101011=\text{000011}.00110001$	3	$0.00110001\ast 000110=\text{000001}.00100110$	1
20	$010100\ast 0.00101011=\text{000011}.01011100$	3	$0.01011100\ast 000110=\text{000010}.00101000$	2
21	$010101\ast 0.00101011=\text{000011}.10000111$	3	$0.10000111\ast 000110=\text{000011}.00101010$	3
22	$010110\ast 0.00101011=\text{000011}.10110010$	3	$0.10110010\ast 000110=\text{000100}.00101100$	4
23	$010111\ast 0.00101011=\text{000011}.11011101$	3	$0.11011101\ast 000110=\text{000101}.00101110$	5
24	$011000\ast 0.00101011=\text{000100}.00001000$	4	$0.00001000\ast 000110=\text{000000}.00110000$	0
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
63	$111111\ast 0.00101011=\text{001010}.10010101$	10	$0.10010101\ast 000110=\text{000011}.01111110$	3

Table 2. Table 2: Systems Tested

Processor	Microarchitecture	Compilers
Intel i7-6700	Skylake (x64)	GCC 6.2; Clang 4.0	default platform
Intel i7-4770	Haswell (x64)	GCC 5.4; Clang 3.8
AMD Ryzen 7 1700X	Zen (x64)	GCC 7.2; Clang 4.0
POWER8	POWER8 Murano	GCC 5.4; Clang 3.8
AMD Opteron A1100	ARM Cortex A57 (Aarch64)	GCC 6.2; Clang 4.0

Table 3. Table 3: Processor cycles (in gigacycles) to determine the number of primes less than 40000, 1000 times, using unsigned 32-bit computations. Branchful and branchless are libdivide alternatives. Note that libdivide was only available for the x64 systems as it uses platform-specific optimizations. POWER8 results are in user CPU seconds. Boldfacing indicates the fastest approach.

Algorithm	Skylake		Haswell		Ryzen		ARM		POWER8
	GCC	Clang	GCC	Clang	GCC	Clang	GCC	Clang	GCC	Clang
division instruction	72	72	107	107	131	131	65	65	18	17
branchful	46	88	56	98	59	71	–	–	–	–
branchless	35	35	36	37	34	37	–	–	–	–
LKK	18	18	18	18	17	18	27	27	8.7	8.0
GM	24	27	27	28	27	32	36	37	10	11
GM/LKK	1.33	1.50	1.50	1.55	1.59	1.77	1.33	1.37	1.15	1.38

Equations

$$n*c=\stackrel{N=6}{\overbrace{111111}}*0.\underset{\_}{\underbrace{\stackrel{N=6}{\overbrace{001010}}\stackrel{L=2}{\overbrace{11}}}}F=N+L=8=\stackrel{\mathrm{quotient}:N\mathrm{bits}}{\overbrace{0001010}}.\stackrel{(c*n)\mathrm{mod}{2}^F}{\overbrace{10010101}}.$$

$$((c*n)\mathrm{mod}{2}^F)*d=\stackrel{(c*n)\mathrm{mod}{2}^F}{\overbrace{10010101}}*\stackrel{d:N\mathrm{bits}}{\overbrace{000110}}=\stackrel{\mathrm{remainder}:N\mathrm{bits}}{\overbrace{0000011}}.\cdots .$$

$$2^{N+L}\le c*d\le 2^{N+L}+2^L$$

$$c*(n\mathrm{mod}d)\le (c*n)\mathrm{mod}{2}^{N+L}\le c*(n\mathrm{mod}d)+2^L(n\mathrm{div}d)<2^{N+L}$$

$$c*r+2^{L}q$$

$$c*r+2^{L}q<2^{N+L}.$$

$$0\le c*d*q-2^{N+L}q\le 2^{L}q.$$

$$c*r\le c*n-2^{N+L}q\le 2^{L}q+c*r$$

$$\frac{1}{d}\le \frac{c}{2^{N+L}}\le \frac{1}{d}+\frac{1/d}{2^N}$$

$$n\mathrm{mod}d=\left(\left((c*n)\mathrm{mod}{2}^{N+L}\right)*d\right)\mathrm{div}{2}^{N+L}$$

$$2^{N-1+L}<c*d<2^{N-1+L}+2^L,$$

$$\mu$$

Full text

Currently, the remainder of the division by a constant is computed from the quotient by a multiplication and a subtraction. But if just the remainder is desired and the quotient is unneeded, this may be suboptimal. We present a generally applicable algorithm to compute the remainder more directly. Specifically, we use the fractional portion of the product of the numerator and the inverse of the divisor. On this basis, we also present a new, simpler divisibility algorithm to detect nonzero remainders.

We also derive new tight bounds on the precision required when representing the inverse of the divisor. Furthermore, we present simple C implementations that beat the optimized code produced by state-of-art C compilers on recent x64 processors (e.g., Intel Skylake and AMD Ryzen), sometimes by more than 25%. On all tested platforms including 64-bit ARM and POWER8, our divisibility-test functions are faster than state-of-the-art Granlund-Montgomery divisibility-test functions, sometimes by more than 50%.

\jnlcitation\cname

, , and , \ctitleFaster Remainder by Direct Computation, \cvol2018.

††articletype: Research Article

1 Introduction

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 1. We find similar latencies in the popular ARM processors 2. 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.

2 Related Work

Jacobsohn3 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, Vowels4 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 Montgomery5 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 Werbrouck6, Warren7 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.8, Granlund and Montgomery5 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 Montgomery5 state that the remainder, if desired, can be computed by an additional multiplication and subtraction. Warren7 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 Jacobsohn3 and Vowels4 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 literature9, 10, 11. Moreover, many researchers12, 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).

3 Computing the Remainder Directly

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

[TABLE]

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:

[TABLE]

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**.**

Given $d\in[1,2^{N})$, and non-negative integers $c,L$ such that

[TABLE]

then

[TABLE]

for all $n\in[0,2^{N})$.

Proof 3.2**.**

We can write $n$ uniquely as $n=q\operatorname{\ast}d+r$ for some integers $q$ and $r$ where $q\geq 0$ and $r\in[0,d)$. We assume that $2^{N+L}\leq c\operatorname{\ast}d\leq 2^{N+L}+2^{L}$.

We begin by showing that $c\operatorname{\ast}r+2^{L}q<2^{N+L}$. Because $c\operatorname{\ast}d\leq 2^{N+L}+2^{L}$, we have that

[TABLE]

Because $n<2^{N}$ and $r<d$, we have that $n+2^{N}r<2^{N}d$ which shows that

[TABLE]

We can rewrite our assumption $2^{N+L}\leq c\operatorname{\ast}d\leq 2^{N+L}+2^{L}$ as $0\leq c\operatorname{\ast}d-2^{N+L}\leq 2^{L}$. Multiplying throughout by the non-negative integer $q$, we get

[TABLE]

After adding $c\operatorname{\ast}r$ throughout, we get

[TABLE]

where we used the fact that $c\operatorname{\ast}d\operatorname{\ast}q+c\operatorname{\ast}r=c\operatorname{\ast}n$. So we have that $c\operatorname{\ast}n-2^{N+L}q\in[c\operatorname{\ast}r,2^{L}q+c\operatorname{\ast}r]$. We already showed (see Equation 1) that $2^{L}q+c\operatorname{\ast}r$ is less than $2^{N+L}$ so that $c\operatorname{\ast}n-2^{N+L}q\in[0,2^{N+L})$. Thus we have that $c\operatorname{\ast}n-2^{N+L}q=(c\operatorname{\ast}n)\bmod 2^{N+L}$ because (in general and by definition) if $p-kQ\in[0,y)$ for some $y\leq Q$, then $p\bmod Q=p-kQ$. Hence, we have that $(c\operatorname{\ast}n)\bmod 2^{N+L}\in[c\operatorname{\ast}r,2^{L}q+c\operatorname{\ast}r]$. This completes the proof.

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.

Theorem 3.3**.**

Given $d\in[1,2^{N})$, and non-negative integers $c,L$ such that

[TABLE]

then

[TABLE]

for all $n\in[0,2^{N})$.

Proof 3.4**.**

We can write $n$ uniquely as $n=q\operatorname{\ast}d+r$ where $q\geq 0$ and $r\in[0,d)$. By Lemma 3.1, we have that $c\operatorname{\ast}n\bmod 2^{N+L}\in[c\operatorname{\ast}r,c\operatorname{\ast}r+2^{L}q]$ for all $n\in[0,2^{N})$.

We want to show that if we multiply any value in $[c\operatorname{\ast}r,c\operatorname{\ast}r+2^{L}q]$ by $d$ and divide it by $2^{N+L}$, then we get $r$. That is, if $y\in[c\operatorname{\ast}r,c\operatorname{\ast}r+2^{L}q]$, then $d\operatorname{\ast}y\in[2^{N+L}r,2^{N+L}(r+1))$. We can check this inclusion using two inequalities:

• •

($d\operatorname{\ast}y\geq 2^{N+L}r$) It is enough to show that $c\operatorname{\ast}d\operatorname{\ast}r\geq 2^{N+L}r$ which follows since $c\operatorname{\ast}d\geq 2^{N+L}$ by one of our assumptions.

• •

($d\operatorname{\ast}y<2^{N+L}(r+1)$) It is enough to show that $d\operatorname{\ast}(c\operatorname{\ast}r+2^{L}q)<2^{N+L}(r+1)$. Using the assumption that $c\operatorname{\ast}d\leq 2^{N+L}+2^{L}$, we have that $d\operatorname{\ast}(c\operatorname{\ast}r+2^{L}q)\leq 2^{N+L}r+2^{L}r+2^{L}d\operatorname{\ast}q=2^{N+L}r+2^{L}n$. Since $n<2^{N}$, we finally have $d\operatorname{\ast}(c\operatorname{\ast}r+2^{L}q)<2^{N+L}(r+1)$ as required.

This concludes the proof.

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_{\_}2(d)$.

Lemma 3.5**.**

Given a divisor $d\in[1,2^{N})$, if we set $c=\left\lceil 2^{N+L}/d\right\rceil$, then

• •

$c\operatorname{\ast}d=2^{N+L}+d-(2^{N+L}\bmod d)$* when $d$ is not a power of two,*

• •

and $c\operatorname{\ast}d=2^{N+L}$ when $d$ is a power of two.

Proof 3.6**.**

The case when $d$ is a power of two follows by inspection, so suppose that $d$ is not a power of two. We seek to round $2^{N+L}$ up to the next multiple of $d$. The previous multiple of $d$ is smaller than $2^{N+L}$ by $2^{N+L}\bmod d$. Thus we need to add $d-2^{N+L}\bmod d$ to $2^{N+L}$ to get the next multiple of $d$.

Example 3.7**.**

Consider Table 1 where we divide by $d=6$ and we want to support the numerators $n$ between 0 and $64=2^{6}$ so that $N=6$. It is enough to pick $F\geq N+\log_{\_}2(d)\approx 8.58$ or $F=9$ but we can do better. According to Algorithm 2, the number of fractional bits $F=N+L=6+L$ must satisfy $d\leq(2^{6+L}\bmod d)+2^{L}$. Picking $L=0$ and $F=6$ does not work since $2^{6}\bmod 6+1=5$. Picking $L=1$ also does not work since $2^{7}\bmod 6+2=4$. Thus we need $L=2$ and $F=8$, at least. So we can pick $c=\left\lceil 2^{8}/6\right\rceil=43$. In binary, representing 43 with 8 fractional bits gives $0.00101011$, as in Table 1. Let us divide $63$ by $6$. The quotient is $(63\operatorname{\ast}43)\operatorname*{\mathrm{div}}2^{8}=10$. The remainder is $(((63\operatorname{\ast}43)\bmod 2^{8})\operatorname{\ast}6)\operatorname*{\mathrm{div}}2^{8}=3$.

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}$).

3.1 Directly Computing Signed Remainders

Having established that we could compute the remainder in the unsigned case directly, without first computing the quotient, we proceed to establish the same in the signed integer case. We assume throughout that the processor represents signed integers in $[-2^{N-1},2^{N-1})$ using the two’s complement notation. We assume that the integer $N\geq 1$.

Though the quotient and the remainder have a unique definition over positive integers, there are several valid ways to define them over signed integers. We adopt a convention regarding signed integer arithmetic that is widespread among modern computer languages, including C99, Java, C#, Swift, Go, and Rust. Following Granlund and Montgomery 5, we let $\operatorname{\mathrm{trunc}}(v)$ be $v$ rounded towards zero: it is $\left\lfloor v\right\rfloor$ if $v\geq 0$ and $\left\lceil v\right\rceil$ otherwise. We use “$\operatorname*{\operatorname*{\mathrm{div}}}$” to denote the signed integer division defined as $n\operatorname*{\operatorname*{\mathrm{div}}}d\equiv\operatorname{\mathrm{trunc}}(n/d)$ and “$\operatorname*{\bmod}$” to denote the signed integer remainder defined by the identity $n\operatorname*{\bmod}d\equiv n-\operatorname{\mathrm{trunc}}(n/d)\operatorname{\ast}d$. Changing the sign of the divisor changes the sign of the quotient but the remainder is insensitive to the sign of the divisor. Changing the sign of the numerator changes the sign of both the quotient and the remainder: $(-n)\operatorname*{\operatorname*{\mathrm{div}}}d=-(n\operatorname*{\operatorname*{\mathrm{div}}}d)$ and $(-n)\operatorname*{\bmod}d=-(n\operatorname*{\bmod}d)$.

Let $\operatorname*{\textrm{lsb}}_{\_}K(n)$ be the function that selects the $K$ least significant bits of an integer $n$, zeroing others. The result is always non-negative (in $[0,2^{K})$) in our work. Whenever $2^{K}$ divides $n$, whether $n$ is positive or not, we have $\operatorname*{\textrm{lsb}}_{\_}K(n)=0=n\bmod 2^{K}$.

Remark 3.8**.**

Suppose $2^{K}$ does not divide $n$, then $n\bmod 2^{K}=2^{K}-\operatorname*{\textrm{lsb}}_{\_}K(n)$ when the integer $n$ is negative, and $n\bmod 2^{K}=\operatorname*{\textrm{lsb}}_{\_}K(n)$ when it is positive. Thus we have $\operatorname*{\textrm{lsb}}_{\_}K(n)+\operatorname*{\textrm{lsb}}_{\_}K(-n)=2^{K}$ whenever $2^{K}$ does not divide $n$.

We establish a few technical results before proving Theorem 3.15. We believe it is novel.

Lemma 3.9**.**

Given $d\in[1,2^{N-1})$, and non-negative integers $c,L$ such that

[TABLE]

we have that $2^{N-1+L}$ cannot divide $c\operatorname{\ast}n$ for any $n\in[-2^{N-1},0)$.

Proof 3.10**.**

First, we prove that $2^{L}$ cannot divide $c$. When $2^{L}$ divides $c$, setting $c=\alpha 2^{L}$ for some integer $\alpha$, we have that $2^{N-1}<\alpha d<2^{N-1}+1$, but there is no integer between $2^{N-1}$ and $2^{N-1}+1$.

Next, suppose that $n\in[-2^{N-1},0)$ and $2^{N-1+L}$ divides $c\operatorname{\ast}n$. Since $2^{N-1+L}$ divides $c\operatorname{\ast}n$, we know that the prime factorization of $c\operatorname{\ast}n$ has at least $N-1+L$ copies of 2. Within the range of $n$ ($[-2^{N-1},0)$) at most $N-1$ copies of 2 can be provided. Obtaining the required $N-1+L$ copies of 2 is only possible when $n=-2^{N-1}$ and $c$ provides the remaining copies—so $2^{L}$ divides $c$. But that is impossible.

Lemma 3.11**.**

Given $d\in[1,2^{N-1})$, and non-negative integers $c,L$ such that

[TABLE]

for all integers $n\in(0,2^{N-1}]$ and letting $y=\operatorname*{\textrm{lsb}}_{\_}{N-1+L}(c\operatorname{\ast}n)$, we have that $(y\operatorname{\ast}d)\bmod 2^{N-1+L}>0$.

Proof 3.12**.**

Write $n=q\operatorname{\ast}d+r$ where $r\in[0,d)$.

Since $c\operatorname{\ast}n$ is positive we have that $y=(c\operatorname{\ast}n)\bmod 2^{N-1+L}=\operatorname*{\textrm{lsb}}_{\_}{N-1+L}(c\operatorname{\ast}n)$. (See Remark 3.8.)

Lemma 3.1 is applicable, replacing $N$ with $N-1$. We have a stronger constraint on $c\operatorname{\ast}d$, but that is not harmful. Thus we have $y\in[c\operatorname{\ast}r,c\operatorname{\ast}r+2^{L}q]$.

We proceed as in the proof of Theorem 3.3. We want to show that $d\operatorname{\ast}y\in(2^{N-1+L}r,2^{N-1+L}(r+1))$.

• •

($d\operatorname{\ast}y>2^{N-1+L}r$) We have $y\geq c\operatorname{\ast}r$. Multiplying throughout by $d$, we get $d\operatorname{\ast}y\geq c\operatorname{\ast}d\operatorname{\ast}r>2^{N-1+L}r$, where we used $c\operatorname{\ast}d>2^{N-1+L}$ in the last inequality.

• •

($d\operatorname{\ast}y<2^{N-1+L}(r+1)$) We have $y\leq c\operatorname{\ast}r+2^{L}q$. Multiplying throughout by $d$, we get $d\operatorname{\ast}y\leq c\operatorname{\ast}d\operatorname{\ast}r+2^{L}q\operatorname{\ast}d<2^{N-1+L}r+2^{L}n$. Because $n\leq 2^{N-1}$, we have the result $d\operatorname{\ast}y<2^{N-1+L}(r+1)$.

Thus we have $d\operatorname{\ast}y\in(2^{N-1+L}r,2^{N-1+L}(r+1))$, which shows that $(y\operatorname{\ast}d)\bmod 2^{N-1+L}>0$. This completes the proof.

Lemma 3.13**.**

Given positive integers $a,b,d$, we have that $\left\lfloor(b-a)\operatorname{\ast}d/b\right\rfloor=d-1-\left\lfloor a\operatorname{\ast}d/b\right\rfloor$ if $b$ does not divide $a\operatorname{\ast}d$.

Proof 3.14**.**

Define $p=a\operatorname{\ast}d/b-\left\lfloor a\operatorname{\ast}d/b\right\rfloor$. We have $p\in(0,1)$. We have that $\left\lfloor(b-a)\operatorname{\ast}d/b\right\rfloor=\left\lfloor d-a\operatorname{\ast}d/b\right\rfloor=\left\lfloor d-\left\lfloor a\operatorname{\ast}d/b\right\rfloor-p\right\rfloor=d-1-\left\lfloor a\operatorname{\ast}d/b\right\rfloor.$

Theorem 3.15**.**

Given $d\in[1,2^{N-1})$, and non-negative integers $c,L$ such that

[TABLE]

let $\mu=\left\lfloor(\operatorname*{\textrm{lsb}}_{\_}{N-1+L}(c\operatorname{\ast}n)\operatorname{\ast}d)/2^{N-1+L}\right\rfloor$ then

• •

$n\operatorname*{\bmod}d=\mu$* for all $n\in[0,2^{N-1})$*

• •

and $n\operatorname*{\bmod}d=\mu-d+1$ for all $n\in[-2^{N-1},0)$.

Proof 3.16**.**

When $n$ is non-negative, then so is $c\operatorname{\ast}n$, and $\operatorname*{\textrm{lsb}}_{\_}{N-1+L}(c\operatorname{\ast}n)$ is equal to $(c\operatorname{\ast}n)\bmod 2^{N-1+L}$; there is no distinction between signed and unsigned $\bmod$, so the result follows by Theorem 3.3, replacing $N$ by $N-1$. (Theorem 3.3 has a weaker constraint on $c\operatorname{\ast}d$.)

Suppose that $n$ is negative ($n\in[-2^{N-1},0)$). By Lemma 3.9, $2^{N-1+L}$ cannot divide $c\operatorname{\ast}n$. Hence, we have that $\operatorname*{\textrm{lsb}}_{\_}{N-1+L}(c\operatorname{\ast}n)=2^{N-1+L}-\operatorname*{\textrm{lsb}}_{\_}{N-1+L}(c\operatorname{\ast}(-n))$ by Remark 3.8. Thus we have

[TABLE]

Hence we have $\mu-d+1=-((-n)\bmod d)=n\bmod d$, which concludes the proof.

We do not need to be concerned with negative divisors since $n\operatorname*{\bmod}d=n\operatorname*{\bmod}-d$ for all integers $n$.

We can pick $c,L$ in a manner similar to the unsigned case. We can choose $c=\left\lfloor 2^{F}/d\right\rfloor+1$ and let $F=N-1+L$ where $L$ is an integer such that $2^{N-1+L}<c\operatorname{\ast}d<2^{N-1+L}+2^{L}$. With this choice of $c$, we have that $c\operatorname{\ast}d=2^{N-1+L}-(2^{N-1+L}\bmod d)+d$. Thus we have the constraint $-(2^{N-1+L}\bmod d)+d<2^{L}$ on $L$. Because $-(2^{N-1+L}\bmod d)+d\in[1,d]$, it suffices to pick $L$ large enough so that $2^{L}>d$. Thus any $L>\log_{\_}2d$ would do, and hence $F>N+\log_{\_}2(d)$ is sufficient. It is not always best to pick $L$ to be minimal: it could be convenient to pick $L=N+1$.

3.2 Fast Divisibility Check with a Single Multiplication

Following earlier work by Artzy et al.8, Granlund and Montgomery5 describe how we can check quickly whether an unsigned integer is divisible by a constant, without computing the remainder. We summarize their approach before providing an alternative. Given an odd divisor $d$, we can find its (unique) multiplicative inverse $\bar{d}$ defined as $d\operatorname{\ast}\bar{d}\bmod 2^{N}=1$. The existence of a multiplicative inverse $\bar{d}$ allows us to quickly divide an integer $n$ by $d$ when it is divisible by $d$, if $d$ is odd. It suffices to multiply $n=a\operatorname{\ast}d$ by $\bar{d}$: $n\operatorname{\ast}\bar{d}\bmod 2^{N}=a\operatorname{\ast}(d\operatorname{\ast}\bar{d})\bmod 2^{N}=a\bmod 2^{N}=n\operatorname*{\mathrm{div}}d$. When the divisor is $2^{K}d$ for $d$ odd and $n$ is divisible by $2^{K}d$, then we can write $n\operatorname*{\mathrm{div}}(2^{K}d)=(n\operatorname*{\mathrm{div}}2^{K})\operatorname{\ast}\bar{d}\bmod 2^{N}$. As pointed out by Granlund and Montgomery, this observation can also enable us to quickly check whether a number is divisible by $d$. If $d$ is odd and $n\in[0,2^{N})$ is divisible by $d$, then $n\operatorname{\ast}\bar{d}\bmod 2^{N}\in[0,\left\lfloor(2^{N}-1)/d\right\rfloor]$. Otherwise $n$ is not divisible by $d$. Thus, when $d$ is odd, we can check whether any integer in $[0,2^{N})$ is divisible by $d$ with a multiplication followed by a comparison. When the divisor is even ($2^{K}d$), then we have to check that $n\operatorname{\ast}\bar{d}\bmod 2^{N}\in[0,2^{K}\operatorname{\ast}\left\lfloor(2^{N}-1)/d\right\rfloor]$ and that $n$ is divisible by $2^{K}$ (i.e., $n\bmod 2^{K}=0$). We can achieve the desired result by computing $n\operatorname{\ast}\bar{d}$, rotating the resulting word by $K$ bits and comparing the result with $\left\lfloor(2^{N}-1)/d\right\rfloor$.

Granlund and Montgomery can check that an unsigned integer is divisible by another using as little as one multiplication and comparison when the divisor is odd, and a few more instructions when the divisor is even. Yet we can always check the divisibility with a single multiplication and a modulo reduction to a power of two—even when the divisor is even because of the following proposition. Moreover, a single precomputed constant ($c$) is required.

Proposition 3.17**.**

Given $d\in[1,2^{N})$, and non-negative integers $c,L$ such that $2^{N+L}\leq c\operatorname{\ast}d\leq 2^{N+L}+2^{L}$ then given some $n\in[0,2^{N})$, we have that $d$ divides $n$ if and only if $(c\operatorname{\ast}n)\bmod 2^{N+L}<c$.

Proof 3.18**.**

We have that $d$ divides $n$ if and only if $n\bmod d=0$. By Lemma 3.1, we have that $c\operatorname{\ast}(n\bmod d)\leq(c\operatorname{\ast}n)\bmod 2^{N+L}\leq c\operatorname{\ast}(n\bmod d)+2^{L}(n\operatorname*{\mathrm{div}}d)$. We want to show that $n\bmod d=0$ is equivalent to $(c\operatorname{\ast}n)\bmod 2^{N+L}<c$.

Suppose that $n\bmod d=0$, then we have that $(c\operatorname{\ast}n)\bmod 2^{N+L}\leq 2^{L}(n\operatorname*{\mathrm{div}}d)$. However, by our constraints on $c$, we have that $c\geq 2^{N+L}/d>2^{L}(n\operatorname*{\mathrm{div}}d)$. Thus, if $n\bmod d=0$ then $(c\operatorname{\ast}n)\bmod 2^{N+L}<c$.

Suppose that $(c\operatorname{\ast}n)\bmod 2^{N+L}<c$, then because $c\operatorname{\ast}(n\bmod d)\leq(c\operatorname{\ast}n)\bmod 2^{N+L}$, we have that $c\operatorname{\ast}(n\bmod d)<c$ which implies that $n\bmod d=0$. This completes the proof.

Thus if we have a reciprocal $c=\left\lceil 2^{F}/d\right\rceil$ with $F=N+L$ large enough to compute the remainder exactly (see Algorithm 2), then $(c\operatorname{\ast}n)\bmod 2^{F}<c$ if and only if $n$ is divisible by $d$. We do not need to pick $F$ as small as possible. In particular, if we set $c=\left\lceil 2^{2N}/d\right\rceil$, then $(c\operatorname{\ast}n)\bmod 2^{2N}<c$ if and only if $n$ is divisible by $d$.

Remark 3.19**.**

We can extend our fast divisibility check to the signed case. Indeed, we have that $d$ divides $n$ if and only if $|d|$ divides $|n|$. Moreover, the absolute value of any $N$-bit negative integer can be represented as an $N$-bit unsigned integer.

4 Software Implementation

Using the C language, we provide our implementations of the 32-bit remainder computation (i.e., a % d) in Figs. LABEL:fig:codemod and 1 for unsigned and signed integers. In both case, the programmer is expected to precompute the constant c. For simplicity, the code shown here explicitly does not handle the divisors $d\in\{-1,0,1,-2^{31}\}$.

For the x64 platforms, we provide the instruction sequences in assembly code produced by GCC and Clang for computing $n\bmod 95$ in Fig. 2; in the third column, we provide the x64 code produced with our approach after constant folding. Our approach generates about half as many instructions.

In Fig. 3, we make the same comparison on the 64-bit ARM platform, putting side-by-side compiler-generated code for the Granlund-Montgomery-Warren approach with code generated from our approach. As a RISC processor, ARM does not handle most large constants in a single machine instruction, but typically assembles them from 16-bit quantities. Since the Granlund-Montgomery-Warren algorithm requires only 32-bit constants, two 16-bit values are sufficient whereas our approach relies on 64-bit quantities and thus needs four 16-bit values. The ARM processor also has a “multiply-subtract” instruction that is particularly convenient for computing the remainder from the quotient. Unlike the case with x64, our approach does not translate into significantly fewer instructions on the ARM platform.

These code fragments show that a code-size saving is achieved by our approach on x64 processors, compared to the approach taken by the compilers. We verify in § 5 that there is also a runtime advantage.

4.1 Divisibility

We are interested in determining quickly whether a 32-bit integer $d$ divides a 32-bit integer $n$—faster than by checking whether the remainder is zero. To the best of our knowledge, no compiler includes such an optimization, though some software libraries provide related fast functions.bbbhttps://gmplib.com We present the code for our approach (LKK) in Fig. 4, and our implementation of the Granlund-Montgomery approach (GM) in Fig. 5.

5 Experiments

Superiority over the Granlund-Montgomery-Warren approach might depend on such CPU characteristics as the relative speeds of instructions for integer division, 32-bit integer multiplication and 64-bit integer division. Therefore, we tested our software on several x64 platforms and on ARMccc With GCC 4.8 on the ARM platform we observed that, for many constant divisors, the compiler chose to generate a udiv instruction instead of using the Granlund-Montgomery code sequence. This is not seen for GCC 6.2. and POWER8 servers, and relevant details are given in Table 2. The choice of multiplication instructions and instruction scheduling can vary by compiler, and thus we tested using various versions of GNU GCC and LLVM’s Clang. For brevity we primarily report results from the Skylake platform, with comments on points where the other platforms were significantly different. For the Granlund-Montgomery-Warren approach with compile-time constants, we use the optimized divide and remainder operations built into GCC and Clang.

We sometimes need to repeatedly divide by a constant that is known only at runtime. In such instances, an optimizing compiler may not be helpful. Instead a programmer might rely on a library offering fast division functions. For runtime constants on x64 processors, we use the libdivide librarydddhttp://libdivide.com as it provides a well-tested and optimized implementation.

On x64 platforms, we use the compiler flags -O3 -march=native; on ARM we use -O3 -march=armv8-a and on POWER8 we use -O3 -mcpu=power8. Some tests have results reported in wall-clock time, whereas in other tests, the Linux perf stat command was used to obtain the total number of processor cycles spent doing an entire benchmark program. To ease reproducibility, we make our benchmarking software and scripts freely available.eeehttps://github.com/lemire/constantdivisionbenchmarks

5.1 Beating the Compiler

We implement a 32-bit linear congruential generator14 that generates random numbers according to the function $X_{\_}{n+1}=(a\operatorname{\ast}X_{\_}n+b)\bmod d$, starting from a given seed $X_{\_}0$. Somewhat arbitrarily, we set the seed to 1234, we use 31 as the multiplier ($a=31$) and the additive constant is set to 27961 ($b=27961$). We call the function 100 million times, thus generating 100 million random numbers. The divisor $d$ is set at compile time. See Fig. 6. In the signed case, we use a negative multiplier ($a=-31$).

Because the divisor is a constant, compilers can optimize the integer division using the Granlund-Montgomery approach. We refer to this scenario as the compiler case. To prevent the compiler from proceeding with such an optimization and force it to repeatedly use the division instruction, we can declare the variable holding the modulus to be volatile (as per the C standard). We refer to this scenario as the division instruction case. In such cases, the compiler is not allowed to assume that the modulus is constant—even though it is. We verified that the assembly generated by the compiler includes the division instruction and does not include expensive operations such as memory barriers or cache flushes. We verified that our wall-clock times are highly repeatablefffFor instance, we repeated tests 20 times for 9 divisors in Figs. 7abcd and 8ab, and we observed maximum differences among the 20 trials of 4.8 %, 0.3 %, 0.7 %, 0.0 %, 0.8 % and 0.9 %, respectively. .

We present our results in Fig. 7 where we compare with our alternative. In all cases, our approach is superior to the code produced by the compiler, except for powers of two in the case of GCC. The benefit of our functions can reach 30%.

The performance of the compiler (labelled as compiler) depends on the divisor for both GCC and Clang, though Clang has greater variance. The performance of our approach is insensitive to the divisor, except when the divisor is a power of two.

We observe that in the unsigned case, Clang optimizes very effectively when the divisor is a small power of two. This remains true even when we disable loop unrolling (using the -fno-unroll-loops compiler flag). By inspecting the produced code, we find that Clang (but not GCC) optimizes away the multiplication entirely in the sense that, for example, $X_{\_}{n+1}=(31\operatorname{\ast}X_{\_}n+27961)\bmod 16$ is transformed into $X_{\_}{n+1}=\operatorname*{\textrm{lsb}}_{\_}{4}(9-X_{\_}n)$. We find it interesting that these optimizations are applied both in the compiler functions as well as in our functions. Continuing with the unsigned case, we find that Clang often produces slightly more efficient compiled code than GCC for our functions, even when the divisor is not a power of two: compare Fig. 7a with Fig. 7c. However, these small differences disappear if we disable loop unrolling.

Yet, GCC seems preferable in the signed benchmark: in Figs. 7b and 7d, Clang is slightly less efficient than GCC, sometimes requiring $0.5\text{\,}\mathrm{s}$ to complete the computation whereas GCC never noticeably exceeds $0.4\text{\,}\mathrm{s}$.

For comparison, Fig. 8 shows how the Ryzen, POWER8 and ARM processors perform on unsigned computations.The speed of the hardware integer-division instruction varies, speeding up at $d=8$ and again at $d=32$ for Ryzen and $d=4$, 16, 64, 256 and 1024 for ARM. The gap between hardware integer division and Granlund-Montgomery (compiler) is less on Ryzen, POWER8 and ARM than on Skylake; for some divisors, there is little benefit to using compiler on POWER8 and ARM. On x64 platforms, our approach continues to be significantly faster than hardware integer division for all divisors.

On ARM, the performance is limited when computing remainders using our approach. Unlike x64 processors, ARM processors require more than one instruction to load a constant such as the reciprocal ($c$), but that is not a concern in this instance since the compiler loads $c$ into a register outside of the loop. We believe that the reduced speed has to do with the performance of the multiplication instructions of our Cortex A57 processor 2. To compute the most significant 64 bits of a 64-bit product as needed by our functions, we must use the multiply-high instructions (umulh and smulh), but they require six cycles of latency and they prevent the execution of other multi-cycle instructions for an additional three cycles. In contrast, multiplication instructions on x64 Skylake processors produce the full 128-bit product in three cycles. Furthermore, our ARM processor has a multiply-and-subtract instruction with a latency of three cycles. Thus it is advantageous to rely on the multiply-and-subtract instruction instead of the multiply-high instruction. Hence, it is faster to compute the remainder from the quotient by multiplying and subtracting ($r=n-(n\operatorname*{\mathrm{div}}d)\operatorname{\ast}d$). Furthermore, our ARM processor has fast division instructions: the ARM optimization manual for Cortex A57 processors indicates that both signed and unsigned division require between 4 and 20 cycles of latency 2 whereas integer division instructions on Skylake processors (idiv and div) have 26 cycles of latency for 32-bit registers 1. Even if we take into account that division instructions on ARM computes solely the quotient, as opposed to both the quotient and remainder on x64, it seems that the ARM platform has a competitive division latency. Empirically, the division instruction on ARM is often within 10% of the Granlund-Montgomery compiler optimization (Fig. 8b) whereas the compiler optimization is consistently more than twice as fast as the division instruction on a Skylake processor (see Fig. 7a).

Results for POWER8 are shown in Figs. 8c and 8d. Our unsigned approach is better than the compiler’s; indeed the compiler would sometimes have done better to generate a divide instruction than use the Granlund-Montgomery approach. For our signed approach, both GCC and Clang had trouble generating efficient code for many divisors.

As with ARM, code generated for POWER8 also deals with 64-bit constants less directly than x64 processors. If not in registers, POWER8 code loads 64-bit constants from memory, using two operations to construct a 32-bit address that is then used with a load instruction. In this benchmark, however, the compiler keeps 64-bit constants in registers. Like ARM, POWER8 has instructions that compute the upper 64 bits of a 64-bit product. The POWER8 microarchitecture15 has good support for integer division: it has two fixed-point pipelines, each containing a multiplier unit and a divider unit. When the multi-cycle divider unit is operating, fixed-point operations can usually be issued to other units in its pipeline. In our benchmark, dependencies between successive division instructions prevent the processor from using more than one divider. Though we have not seen published data on the actual latency and throughput of division and multiplication on this processor, we did not observe the divisor affecting the division instruction’s speed, at least within the range of 3 to 4096.

Our results suggest that the gap between multiplication and division performance on the POWER8 lies between that of ARM and Intel; the fact that our approach (using 64-bit multiplications) outperforms the compiler’s approach (using 32-bit multiplications) seems to indicate that, unlike ARM, the instruction to compute the most significant bits of a 64-bit product is not much slower than the instruction to compute a 32-bit product.

Looking at Fig. 9, we see how the approaches compare for larger divisors. The division instruction is sometimes the fastest approach on ARM, and sometimes it can be faster than the compiler approach on Ryzen. Overall, our approach is preferred on Ryzen (as well as Skylake and POWER8), but not on ARM.

5.2 Beating the libdivide Library

There are instances when the divisor might not be known at compile time. In such instances, we might use a library such as a libdivide. We once again use our benchmark based on a linear congruential generator using the algorithms, but this time, we provide the divisor as a program parameter.

The libdivide library does not have functions to compute the remainder, so we use its functions to compute the quotient. It has two types of functions: regular "branchful" ones, those that include some branches that depend on the divisor, and branchless ones. In this benchmark, the invariant divisor makes the branches perfectly predictable, and thus the libdivide branchless functions were always slower. Consequently we omit the branchless results.

We present our results in Fig. 10. The performance levels of our functionsggg When 20 test runs were made for 9 divisors, timing results among the 20 never differed by more than 1%. are insensitive to the divisor, and our performance levels are always superior to those of the libdivide functions (by about 15%), except for powers of two in the unsigned case. In these cases, libdivide is faster, but this is explained by a fast conditional code path for powers of two.

5.3 Competing for Divisibility

We adapted a prime-counting benchmark distributed with libdivide, specialized to 32-bit operands. The code determines the number of primes in $[2,40000)$ using a simplistic approach: odd numbers in this range are checked for divisibility by any smaller number that has already been determined to be prime. See Fig. 11. When a number is identified as a prime, we compute its scaled approximate reciprocal ($c$) value, which is repeatedly used in future iterations. In this manner, the computation of $c$ is only done once per prime, and not once per trial division by the prime. A major difference from the benchmark using the linear-congruential generator is that we cycle rapidly between different divisors, making it much more difficult to predict branches in the libdivide functions.

In these tests, we compare libdivide against LKK and GM, the fast divisibility tests whose implementations are shown in § 4.1; see Fig. 4 for LKK and Fig. 5 for GM. Divisibility of a candidate prime is checked either using

• •

libdivide to divide, followed by multiplication and subtraction to determine whether the remainder is nonzero;

• •

the Granlund-Montgomery (GM) divisibility check, as in Fig. 5;

• •

the C % operation, which uses a division instruction;

• •

our LKK divisibility check (Fig. 4).

LKK stores 64 bits for each prime; GM requires an additional 5-bit rotation amount. The division-instruction version of the benchmark only needs to store 32 bits per prime. The libdivide approach requires 72 bits per prime, because we explicitly store the primes.

Instruction counts and execution speed are both important. All else being equal, we would prefer that compilers emit smaller instruction sequences. Using a hardware integer division will yield the smallest code, but this might give up too much speed. In the unsigned case, our LKK has a significant code-size advantage over GM—approximately 3 arithmetic instructions to compute our $c$ versus about 11 to compute their required constant. Both fast approaches use a multiplication and comparison for each subsequent divisibility check. GM requires an additional instruction to rotate the result of the multiplication.

Performance results for the unsigned case are shown in Table 3, showing the total number of processor cycles on each platform from 1000 repetitions of the benchmark. On Skylake, 20 repeated computations yielded cycle-count results within 0.3% of each other. For ARM, results were always within 4%. Initially, Ryzen results would sometimes differ by up to 10% within 20 attempts, even after we attempted to control such factors as dynamic frequency scaling. Thus, rather than reporting the first measurement for each benchmark, the given Ryzen results are the average of 11 consecutive attempts (the basic benchmark was essentially executed 11 000 times). Our POWER8 results (except one outlier) were within 7% of one another over multiple trials and so we averaged several attempts (3 for GCC and 7 for Clang) to obtain each data point. Due to platform constraints, POWER8 results are user-CPU times that matched the wall-clock times.

LKK has a clear speed advantage in all cases, including the POWER8 and ARM platforms. LKK is between 15% to 80% faster than GM. Both GM and LKK always are much faster than using an integer division instruction (up to $7\times$ for Ryzen) and they also outperform the best algorithm in libdivide.

6 Conclusion

To our knowledge, we present the first general-purpose algorithms to compute the remainder of the division by unsigned or signed constant divisors directly, using the fractional portion of the product of the numerator with the approximate reciprocal3, 4. On popular x64 processors (and to a lesser extent on POWER), we can produce code for the remainder of the integer division that is faster than the code produced by well regarded compilers (GNU GCC and LLVM’s Clang) when the divisor constant is known at compile time, using small C functions. Our functions are up to 30% faster and with only half the number of instructions for most divisors. Similarly, when the divisor is reused, but is not a compile-time constant, we can surpass a popular library (libdivide) by about 15% for most divisors.

We can also speed up a test for divisibility. Our approach (LKK) is several times faster than the code produced by popular compilers. It is faster than the Granlund-Montgomery divisibility check5, sometimes nearly twice as fast. It is advantageous on all tested platforms (x64, POWER8 and 64-bit ARM).

Though compilers already produce efficient code, we show that additional gains are possible. As future work, we could compare against more compilers and other libraries. Moreover, various additional optimizations are possible, such as for division by powers of two.

Acknowledgments

The work is supported by the Natural Sciences and Engineering Research Council of Canada under grant RGPIN-2017-03910. The authors are grateful to IBM’s Centre for Advanced Studies — Atlantic and Kenneth Kent for access to the POWER 8 system.