Optimizing integer division

Division is the operation CPUs are worst at. A multiply usually retires in a handful of cycles, often pipelined so a new one can start every cycle, while an integer divide can take several times as long and, on plenty of processors old and new, blocks the pipeline until it's done. Whenever the divisor is fixed — a compile-time constant like x / 10, or at least invariant across many iterations of a loop — that division can be replaced with a multiplication and a shift. This article derives exactly why that works, builds the recipe for picking the "magic" multiplier and shift amount from scratch, and shows how far you can push it in plain C without ever touching assembly.

Division as a Scaled Reciprocal

The one division everybody already replaces with a shift is division by a power of two: n / 8 is n >> 3, exactly, for any unsigned n. The reason is that \(1/8\) has a finite binary expansion. For any other divisor \(d\), \(1/d\) is not exactly representable in binary with a bounded number of digits, but we can still get arbitrarily close: pick a shift amount \(s\) and round \(2^s/d\) up to the nearest integer,

\[m = \left\lceil \frac{2^s}{d} \right\rceil.\]

The claim is that \(\lfloor n\cdot m / 2^s\rfloor\) equals \(\lfloor n/d\rfloor\) for every \(n\) up to some bound that depends on how good the rounding of \(m\) turned out to be. To see why, write \(n = q\cdot d + r\) with \(q=\lfloor n/d\rfloor\) and \(0\le r<d\), and let \(k = m\cdot d - 2^s\) be the small "overshoot" introduced by rounding \(m\) up (\(0\le k<d\), and \(k=0\) exactly when \(d\) divides \(2^s\)). Substituting \(m = (2^s+k)/d\),

\[\frac{n\cdot m}{2^s} = \frac{n}{d} + \frac{n\cdot k}{d\cdot 2^s} = q + \frac{r}{d} + \frac{n\cdot k}{d\cdot 2^s} = q + \frac{r\cdot 2^s + n\cdot k}{d\cdot 2^s}.\]

The floor of the left side equals \(q\) precisely when the fraction on the right stays below 1, i.e. when \(r\cdot 2^s + n\cdot k < d\cdot 2^s\). The worst case is the largest possible remainder, \(r=d-1\), which tightens this to

\[n < \frac{2^s}{k}.\]

That's the whole theorem: choose any shift \(s\), round \(m=\lceil 2^s/d\rceil\), compute \(k=m\cdot d-2^s\), and the multiply-shift trick is exact for every \(n\) below \(2^s/k\). Rounding \(m\) up rather than down is what keeps \(k\) small (at most \(d-1\)), which is why the safe range comes out close to \(2^s/d\) instead of something much smaller.

Picking Good Magic Numbers

With the bound \(n < 2^s/k\) in hand, picking a multiplier is a two-line search: try a shift \(s\), round up to get \(m\), check how much headroom \(k\) leaves, and stop once the resulting range comfortably covers whatever values you actually expect. In code:

uint32_t magic_multiplier(uint32_t d, int s) {
    return (uint32_t)((((uint64_t)1 << s) + d - 1) / d);   /* ceil(2^s / d) */
}

Working through a few common divisors this way gives multipliers that turn up in a lot of embedded code and DSP libraries that have no divide instruction at all:

and each one turns into a single line:

static inline uint32_t div_by_100(uint32_t n) {   /* correct for n < 43690 */
    return (uint32_t)(((uint64_t)n * 5243u) >> 19);
}

Note the widened uint64_t intermediate: even though \(n\) and the multiplier both fit in 32 bits, their product generally doesn't, so the multiplication has to happen at double width before the shift throws away the low bits. This matches the well-known rule of thumb for finding these constants by hand — pick a power of two that, once divided by the divisor, is comfortably larger than the biggest number you need to support — which is really just a conservative version of the exact bound \(n<2^s/k\) derived above.

Scaling to the Full Range

All of the multipliers above were tuned for a specific, bounded range of \(n\). What if you want a divide-by-10 that's correct for the entire 32-bit range, not just up to 262144? The derivation already tells us how: we need \(n<2^s/k\) to cover all of \([0,2^{32})\). Setting \(\ell=\lceil\log_2 d\rceil\) and \(s=32+\ell\), the overshoot satisfies \(k<d\le 2^\ell\), so \(2^s/k > 2^s/2^\ell = 2^{32}\) automatically. In other words, choosing \(s=32+\ell\) always yields a multiplier that works for every possible 32-bit \(n\), for any divisor whatsoever — this is exactly the recipe compilers use to fold a constant division into a multiply.

The catch is that this \(m\) doesn't always fit in 32 bits. Take \(d=7\): \(\ell=3\), \(s=35\), and \(m=\lceil 2^{35}/7\rceil=4908534053\), which needs 33 bits. Assembly-era implementations solve this with a fairly intricate correction step (multiply by \(m-2^{32}\) instead, then add \(n\) back in before the final shift) to stay within a machine word. In C we can sidestep that entirely by simply reaching for a wider integer type. On any compiler that supports GCC's 128-bit extension, __uint128_t comfortably holds the product of a 32-bit \(n\) and a 33-bit \(m\), so the single clean formula keeps working unmodified:

static inline uint32_t div_by_7(uint32_t n) {   /* correct for every 32-bit n */
    return (uint32_t)(((__uint128_t)n * 4908534053ull) >> 35);
}

The bit-fitting problem that the classic algorithm works so hard to avoid is a consequence of doing this in a fixed-width register; it isn't inherent to the math. Widening one step further than strictly necessary removes the need for the correction entirely. It doesn't scale forever — a full 64-bit divide can need a multiplier past 65 bits combined with a 64-bit \(n\), which no longer fits even in 128 bits for the hardest divisors, and at that point the add-back correction from the literature is the more practical route.

Why Your Compiler Already Does This — and When You Still Have to

If you write n / 10 in C with a modern compiler, this is already what gets emitted; no explicit division instruction shows up in the binary. The technique and its careful edge cases (signed divisors, divisors that are one off from a power of two, exact division, and so on) were worked out in full generality by Torbjörn Granlund and Peter Montgomery in Division by Invariant Integers using Multiplication [GranlundMontgomery1994], and GCC has used a version of it since the mid-1990s.

That leaves two situations where doing it yourself still pays off. First, run-time invariant divisors: if a divisor is only known once a function starts (say, it's an argument, not a literal), the compiler generally can't fold it, but if you divide by that same value thousands of times afterwards, computing \(m\) and \(s\) once and reusing them is a legitimate manual optimization. Second, targets without a hardware divide at all: older fixed-point DSPs, small microcontrollers, and plenty of soft-core designs simply have no division instruction to fall back on, so a library or hand-written routine has to implement the multiply-shift sequence explicitly rather than relying on the compiler having hardware to target.

Exact Division by an Odd Constant

There's a completely different trick available when you already know the division has no remainder — for instance, subtracting two pointers into an array of fixed-size structs, where the byte difference is always a multiple of the struct size. For odd \(d\), there's a unique \(d^{-1}\) with \(d\cdot d^{-1}\equiv 1 \pmod{2^N}\) on an \(N\)-bit machine, and multiplying by it recovers the quotient directly:

\[n = q\cdot d \implies q\cdot d\cdot d^{-1} \equiv q \pmod{2^N} \implies n\cdot d^{-1} \bmod 2^N = q,\]

with no shift and no rounding, since \(q\) fits in \(N\) bits and the modular product lands on it exactly. The inverse itself is easy to compute without the extended Euclidean algorithm: start from the fact that any odd \(d\) is its own inverse modulo \(2^3=8\) (check the four odd residues 1, 3, 5, 7 — each squares to 1 mod 8), then apply Newton's iteration

\[x \leftarrow x\cdot(2 - d\cdot x) \pmod{2^N},\]

which doubles the number of correct bits every round, so four rounds take the 3 correct bits at the start past 32:

uint32_t modinv_odd(uint32_t d) {           /* d must be odd */
    uint32_t x = d;                            /* correct mod 2^3 */
    for (int i = 0; i < 4; i++)
        x *= 2 - d * x;                          /* mod 2^32, via natural wraparound */
    return x;                                     /* x * d == 1 (mod 2^32) */
}

uint32_t exact_div(uint32_t n, uint32_t d) {    /* only valid when d divides n exactly */
    int e = 0;
    while ((d & 1) == 0) { d >>= 1; n >>= 1; e++; }
    return n * modinv_odd(d);
}

It's worth being precise about why this only works for exact division, since it's tempting to assume it computes a general quotient. Write \(n = q\cdot d + r\) again. Then

\[n\cdot d^{-1} \bmod 2^N = q + r\cdot d^{-1} \bmod 2^N.\]

When \(r=0\) this is just \(q\). But when \(r\ne 0\), \(r\cdot d^{-1}\bmod 2^N\) is not a small correction term — \(d^{-1}\) is some large, essentially unstructured odd number, so \(r\cdot d^{-1}\bmod 2^N\) lands somewhere unpredictable across the full range, and adding it to \(q\) produces a value with no relationship to \(n/d\). Working modulo 256 with \(d=5\) (whose inverse is 205, since \(5\cdot 205=1025\equiv 1\bmod 256\)) makes this concrete:

Only the exact multiples of 5 land on the true quotient; everything else is scattered across the whole 8-bit range with no useful pattern. This exact-division trick is also the basis of a fast divisibility test: since \(n\cdot d^{-1}\bmod 2^N\) equals the true quotient exactly when \(d\mid n\), and is otherwise some out-of-range value, checking whether the result stays below \(2^N/d\) tells you whether the remainder was zero, without ever computing a remainder.

A Related Idea: Barrett Reduction

The same reciprocal-approximation trick generalizes from division to modular reduction. Barrett reduction computes \(n\bmod p\) by first estimating \(q\approx n/p\) with exactly the multiply-shift scheme above, then recovering the remainder as \(n - q\cdot p\) and correcting by at most one extra subtraction if the estimate landed a touch too low. It shows up constantly in cryptographic and hashing code that needs a fixed modulus applied millions of times, for the same underlying reason: multiplying is cheap, and dividing, even to get a remainder, is not.

One more small application of the same reciprocal idea, for a change of base rather than a modulus: since \(1233/4096\approx\log_{10}2\), the number of decimal digits needed to print an \(n\)-bit binary number can be estimated with (n * 1233) >> 12 — the same shift-a-scaled-reciprocal move, just approximating an irrational constant instead of \(1/d\), and a handy way to size a buffer before doing exactly the kind of radix conversion this whole approach was built to accelerate in the first place.

References

• GranlundMontgomery1994 T. Granlund and P. L. Montgomery (1994) Division by Invariant Integers using Multiplication. Proceedings of PLDI '94. https://gmplib.org/~tege/divcnst-pldi94.pdf