Fast Division With Magic Numbers

optimization article zig2026-04-03

In the rank-select article, the select query walks down a hierarchy of sampled positions, and at every level it divides the position by that level's stride. The strides are picked at initialization, so the compiler sees a runtime divisor and emits a hardware div.

Integer division is considerably more expensive than multiplication on most modern x86-64 and ARM cores. A dependent 64-bit multiply usually takes a few cycles, while 64-bit division has higher latency and much lower throughput. The exact costs vary by microarchitecture and, on some processors, by operand values.

The compiler can avoid the divider when the divisor is known at compile time. Here are the two cases side by side:

export fn divideByTen(n: u64) u64 {
    return n / 10;
}

export fn divideRuntime(n: u64, d: u64) u64 {
    return n / d;
}

Zig (LLVM) compiles divideByTen to mulx and a shift on x86-64, or umulh and a shift on ARM64, with the reciprocal of 10 baked into the function. divideRuntime has to emit div or udiv because d is not known until the function is called.

When the same runtime divisor is reused millions of times, we can compute its reciprocal in runtime once and use this strategy ourselves. Daniel Lemire describes the construction in Faster remainders when the divisor is a constant.

The Idea

Dividing by $d$ is multiplying by $1/d$ . We store $1/d$ as a 64-bit fixed-point fraction

m = \left\lceil \frac{2^{64}}{d} \right\rceil

and then for any $n$ :

\left\lfloor \frac{n}{d} \right\rfloor = \left\lfloor \frac{n \times m}{2^{64}} \right\rfloor

so the division becomes a multiply and taking the high 64 bits of the product.

The required fraction is twice as wide as the integers being divided. A 64-bit fraction covers every u32 numerator and divisor, while u64 needs a 128-bit fraction.

Proof

The same argument works for an $N$ -bit fraction whenever $n, d < 2^{N/2}$ . To keep the notation concrete, the proof below uses the u32 case with $N = 64$ . Lemire, Kaser and Kurz give the general proof in Faster Remainder by Direct Computation.

Let

m = \left\lceil \frac{2^{64}}{d} \right\rceil = \frac{2^{64}}{d} + e, \qquad 0 \le e < 1

The value computed before taking the floor is therefore

\frac{n \times m}{2^{64}} = \frac{n}{d} + \frac{n \times e}{2^{64}}

Call the second term $E$ . Since $e < 1$ ,

0 \le E < \frac{n}{2^{64}}

Both $n$ and $d$ are smaller than $2^{32}$ , so $n \times d < 2^{64}$ . Dividing both sides by $d \times 2^{64}$ gives

\frac{n}{2^{64}} < \frac{1}{d}

and therefore $E < 1/d$ .

Now apply integer division to $n$ :

n = qd + r, \qquad 0 \le r < d

This gives $n/d = q + r/d$ . If $r = 0$ , then $n/d = q$ . Since $E < 1/d \le 1$ ,

q \le \frac{n}{d} + E < q + 1

If $r > 0$ , the distance from $n/d$ to the next integer is

(q + 1) - \left(q + \frac{r}{d}\right) = \frac{d-r}{d}

Here $d-r$ is a positive integer, so the distance is at least $1/d$ . Because $E < 1/d$ , the error cannot reach the next integer. In both cases,

\left\lfloor \frac{n \times m}{2^{64}} \right\rfloor = q = \left\lfloor \frac{n}{d} \right\rfloor

pub fn magic(d: u32) u64 {
    std.debug.assert(d > 1);
    return ~@as(u64, 0) / d + 1;
}

When $d$ doesn't divide $2^{64}$ , $\lfloor (2^{64}-1)/d \rfloor = \lfloor 2^{64}/d \rfloor$ and the $+1$ turns the floor into the ceiling. When $d$ is a power of two, $\lfloor (2^{64}-1)/d \rfloor = 2^{64}/d - 1$ and the $+1$ lands exactly on the integer result, so both cases end at $\lceil 2^{64}/d \rceil$ .

The division itself is one widening multiply:

pub fn fastdiv(n: u32, m: u64) u32 {
    return @intCast((@as(u128, n) * m) >> 64);
}

Modulo

Let's write the division formula again, we're looking for $r$ in this case

n = qd + r, \qquad 0 \le r < d.

Multiplying by the magic reciprocal gives

n \times m = q \times 2^{64} + r \times \frac{2^{64}}{d} + e, \qquad m = \left\lceil \frac{2^{64}}{d} \right\rceil

The high 64 bits of $n \times m$ result give $q$ , which is what fastdiv returns. The low 64 bits are the part left after removing q * 2^64, so they contain the remainder $r$ scaled by roughly $2^{64} / d$ plus some error. To turn that remainder back into $r$ we multiply it by $d$ and divide by $2^{64}$ :

pub fn fastmod(n: u32, m: u64, d: u32) u32 {
    const frac = m *% @as(u64, n);
    return @intCast((@as(u128, frac) * d) >> 64);
}

m *% n is a wrapping multiply in u64, so we'll get only the low 64 bits of the product.

The only problem is that m is rounded up, so the low bits also contain a small rounding error. For u32 inputs that error is too small to change the final floor after multiplying back by d.

Divisibility

Let's write the division formula again, but now we're only checking if $r$ is zero:

n = qd + r, \qquad 0 \le r < d.

After multiplying by the magic reciprocal, the low 64 bits contain the remainder $r$ scaled by roughly $2^{64}/d$ plus some error. The value $m$ is the same scale applied to one remainder step: $\lceil 2^{64}/d \rceil$ .

So comparing the low bits with $m$ is the scaled version of asking whether

\frac{r}{d} < \frac{1}{d}.

Since $r$ is an integer remainder, the only value below the first non-zero remainder is $r = 0$ .

pub fn isDivisible(n: u32, m: u64) bool {
    return m *% @as(u64, n) < m;
}

m *% n is the same wrapping multiply as in fastmod, it gives us only the low 64 bits of the product.

If the code only needs the divisibility check and never the remainder itself. This is cheaper than computing n % d (fastmod) and comparing it to zero.

u64

The same construction one level up: a 128-bit magic number $\lceil 2^{128}/d \rceil$ and the high 64 bits of the product:

pub fn magicU64(d: u64) u128 {
    return ~@as(u128, 0) / d + 1;
}

pub fn fastdivU64(n: u64, m: u128) u64 {
    return @intCast((@as(u256, n) * m) >> 128);
}

pub fn fastmodU64(n: u64, m: u128, d: u64) u64 {
    const frac = m *% @as(u128, n);
    return @intCast((@as(u256, frac) * d) >> 128);
}

pub fn isDivisibleU64(n: u64, m: u128) bool {
    return m *% @as(u128, n) < m;
}

LLVM magic

While testing this function I found out writing the function like this works a little bit faster on my ARM machine:

fn fastdivU64(n: u64, m: u128) u64 {
    const lo = ((m & ~@as(u64, 0)) * @as(u128, n)) >> 64;
    const hi = (m >> 64) * @as(u128, n);
    return @intCast((lo + hi) >> 64);
}

LLVM lowers both versions to the same operations but in a slightly different order for macos target.

export fn fastdivU64_split(n: u64, m: u128) u64 {
    const lo = ((m & ~@as(u64, 0)) * @as(u128, n)) >> 64;
    const hi = (m >> 64) * @as(u128, n);
    return @intCast((lo + hi) >> 64);
}

export fn fastdivU64_native(n: u64, m: u128) u64 {
    return @intCast((@as(u256, n) * m) >> 128);
}

The difference is the order of the two umulh instructions. On my machine the split variant starts with the multiply whose result is needed by cmn:

umulh   x8, x1, x0
umulh   x9, x2, x0
mul     x10, x2, x0
cmn     x8, x10
cinc    x0, x9, hs

The u256 form starts the other one:

umulh   x8, x2, x0
umulh   x9, x1, x0
mul     x10, x2, x0
cmn     x10, x9
cinc    x0, x8, hs

The umulh resut feeding cmn has to be ready at the same time as the mul result. But in u256 version that umulh instruction goes second for some reason. It seems like the scheduler has a better chance of making the right choice when instructions are already in order.

So the generated code has the same number of instructions and does the same thing, but due to all the modern CPU internal optimizations the split version wins.

For x86-64 and linux ARM target there's no difference.

Benchmark

In the loop below every call feeds its result into the next input, so each division waits for the one before it:

const std = @import("std");

fn magic(d: u32) u64 {
    return ~@as(u64, 0) / d + 1;
}

fn fastdiv(n: u32, m: u64) u32 {
    return @intCast((@as(u128, n) * m) >> 64);
}

fn magicU64(d: u64) u128 {
    return ~@as(u128, 0) / d + 1;
}

fn fastdivU64(n: u64, m: u128) u64 {
    return @intCast((@as(u256, n) * m) >> 128);
}

fn fastdivU64Split(n: u64, m: u128) u64 {
    const lo = ((m & ~@as(u64, 0)) * @as(u128, n)) >> 64;
    const hi = (m >> 64) * @as(u128, n);
    return @intCast((lo + hi) >> 64);
}

const rounds = 50_000_000;

fn bench(name: []const u8, comptime f: anytype, args: anytype) !void {
    const big_prime_1 = 0xfcbdea48a9ee3875;
    const big_prime_2 = 0xe706823742bc7ebd;
    
    var x: u64 = big_prime_1;
    var timer = try std.time.Timer.start();
    
    for (0..rounds) |_| {
        x = @call(.auto, f, .{x} ++ args) +% big_prime_2;
    }
    std.mem.doNotOptimizeAway(x);
    
    const ns = @as(f64, @floatFromInt(timer.read()));
    std.debug.print("{s:<20} {d:.2} ns/op\n", .{ name, ns / rounds });
}

pub fn main() !void {
    var seed: u64 = 0xBAD_DEAD_C0DE;
    // Prevent compile time optimizations
    // Without this LLVM will be able to track down this constant down to the bench code
    std.mem.doNotOptimizeAway(&seed);
    var r = std.Random.DefaultPrng.init(seed);
    const random = r.random();

    const d32 = random.intRangeAtMost(u32, 2, 1 << 16);
    const d64 = random.intRangeAtMost(u64, 2, 1 << 16);
    const m32 = magic(d32);
    const m128 = magicU64(d64);

    try bench("u32 div", hwDiv32, .{d32});
    try bench("u32 fastdiv", fd32, .{m32});
    try bench("u64 div", hwDiv64, .{d64});
    try bench("u64 fastdiv u256", fd64, .{m128});
    try bench("u64 fastdiv split", fd64Split, .{m128});
    try bench("u32 fastdiv + magic", fd32WithMagic, .{&d32});
    try bench("u64 fastdiv + magic", fd64WithMagic, .{&d64});
}


fn hwDiv32(x: u64, d: u32) u64 {
    return @as(u32, @truncate(x)) / d;
}

fn fd32(x: u64, m: u64) u64 {
    return fastdiv(@truncate(x), m);
}

fn fd32WithMagic(x: u64, d: *const u32) u64 {
    const divisor = @as(*volatile const u32, @volatileCast(d)).*; // Prevent optimizations
    return fastdiv(@truncate(x), magic(divisor));
}

fn hwDiv64(x: u64, d: u64) u64 {
    return x / d;
}

fn fd64(x: u64, m: u128) u64 {
    return fastdivU64(x, m);
}

fn fd64Split(x: u64, m: u128) u64 {
    return fastdivU64Split(x, m);
}

fn fd64WithMagic(x: u64, d: *const u64) u64 {
    const divisor = @as(*volatile const u64, @volatileCast(d)).*; // Prevent optimizations
    return fastdivU64Split(x, magicU64(divisor));
}

On my ARM Apple M1 I get the following results:

u32 div              3.07 ns/op
u32 fastdiv          1.58 ns/op
u32 fastdiv + magic  1.61 ns/op
u64 div              3.18 ns/op
u64 fastdiv u256     1.92 ns/op
u64 fastdiv split    1.79 ns/op
u64 fastdiv + magic  3.89 ns/op

On my x86-64 Intel i7-6700HQ (Skylake, 2.6 GHz):

u32 div              7.88 ns/op
u32 fastdiv          1.66 ns/op
u32 fastdiv + magic  7.91 ns/op
u64 div              10.62 ns/op
u64 fastdiv u256     2.03 ns/op
u64 fastdiv split    2.04 ns/op
u64 fastdiv + magic  35.84 ns/op

Sandbox numbers move around between runs since the box is shared, while the ratios hold. How much the replacement saves depends on the divider: Apple silicon has an unusually fast one and the multiply still cuts the time in half, and on the sandbox machine a dependent u64 division costs 12 ns against under 2 ns for the multiply.

The + magic rows show where the setup cost matters. For u32 on this ARM machine, building the magic number is still cheap enough that the result stays close to plain fastdiv. That is not true in general, and it is especially not true for u64: the magic number is 128 bits wide, so constructing it inside the hot loop is much slower than just using div. The intended use is still reuse: compute m once when the divisor is chosen, store it next to that divisor, and then every later division is only the multiply path.