I recently found this code snippet in the source code of pic-scale-safe:

fn div_round_by_1023(v: u32) -> u32 {
    let round = 1 << 9;
    let w = v + round;
    ((w >> 10) + w) >> 10
}

This function computes $\mathrm{round}(v/1023)$ for inputs $v<2^{20}+2^9-1$ with just a few bit shifts and additions in 32-bit arithmetic. Quite efficient.

But that's not all. This trick only 21 bits for the intermediate results. Other approaches like the multiply-add method require 31 bits for intermediate results to perform the same rounded division by 1023 in the range $v<2^{20}+2^9-1$. While not the case for division by 1023, in general this trick needs at most one additional bit, which can be the difference between being able to use 32-bit arithmetic or having to resort to 64-bit arithmetic. This is especially important in SIMD code and code the compiler is supposed to auto-vectorize.

As I hinted, this trick doesn't work just for division by 1023. In general, it works for any divisor of the form $2^n-1$ for inputs $v<2^{2n}+2^{n-1}-1$. Here is the generalized version:

fn div_round_by_2pn_m1(v: u32, n: u32) -> u32 {
    let round = 1 << (n - 1);
    let w = v + round;
    ((w >> n) + w) >> n
}

This function returns exactly $\mathrm{round}(v/(2^n-1))$ for all inputs $v<2^{2n}+2^{n-1}-1$. For larger inputs, the results are typically close, but not exact. This makes it an approximation for rounded division by $2^n-1$.

Unfortunately, it's not obvious at all why this approximation is exact for $v<2^{2n}+2^{n-1}-1$, and why it stops working at $v=2^{2n}+2^{n-1}-1$.

In this article, I will answer both questions and generalize the trick to (1) increase the input range and (2) support floor and ceil division as well.

Results

Before I start with the derivations and proofs, here is a summary of the results.

1. $\mathrm{round}(\frac{v}{2^n-1})$ can be approximated using:

   $$\begin{array}{rl}{R}_1& :=\lfloor \frac{v+2^{n-1}}{2^n}\rfloor \\ R_{i+1}& :=\lfloor \frac{R_i+v+2^{n-1}}{2^n}\rfloor \end{array}$$

   The approximation $R_i$ is exact for all inputs $v<2^{in}+2^{n-1}-1$.

2. $\lfloor \frac{v}{2^n-1}\rfloor$ can be approximated using:

   $$\begin{array}{rl}{F}_1& :=\lfloor \frac{v+1}{2^n}\rfloor \\ F_{i+1}& :=\lfloor \frac{F_i+v+1}{2^n}\rfloor \end{array}$$

   The approximation $F_i$ is exact for all inputs $v<2^{in}+2^n-2$.

3. $\lceil \frac{v}{2^n-1}\rceil$ can be approximated using:

   $$\begin{array}{rl}{C}_1& :=\lfloor \frac{v+2^n-1}{2^n}\rfloor \\ C_{i+1}& :=\lfloor \frac{C_i+v+2^n-1}{2^n}\rfloor \end{array}$$

   The approximation $C_i$ is exact for all inputs $v<2^{in}$.

Code Generation

I also implemented a little code gen tool that uses these results to generate Rust code. You can set $n$, the iteration count, the rounding mode, and the integer type all operations will be performed in. Everything (except for the integer type) can also be made variable at runtime by checking the Parameter box.

n:ParameterIteration countIterations:ParameterRounding mode:FloorRoundCeilParameterInteger type:u8u16u32u64u128

/// Returns `round(v / (2^n - 1))`.
///
/// Returned values are correct if both:
/// 1. the approximation is exact: v < 2^(2n) + 2^(n-1) - 1
/// 2. no overflow occurs:         v < 2^32 - 2^(n-1) - round(v/(2^n-1)).
fn div_round_by_2pn_m1(v: u32, n: u32) -> u32 {
    debug_assert!(n != 0, "Division by zero");
    let round = 1 << (n - 1);
    let w = v + round;
    ((w >> n) + w) >> n
}

Proving correctness

Formally, the approximation is defined as follows: Let $v\in \mathbb{N},n\in {\mathbb{N}}_1$, then:

To capture how good the approximation is, I will introduce an error term $\delta$ defined as follows:

Consequently, the approximation is exact for an input $v$ iff $\delta =0$.

Now, I will write $v$ a bit differently. Let $v=a(2^n-1)+b+2^{n-1}$ for $a\in \mathbb{Z},b\in {\mathbb{N}}_0,b<2^n-1$. While a bit unusual, it should be obvious that all integers can be uniquely expressed in this form. I choose this form because it has the nice property that $\mathrm{round}(v/(2^n-1))=a+1$.

This makes it possible to simplify the error term:

Now it is easy to show that $-2^n\le b-a<2^n⟹\delta =0$.

All that is left to do, is to show that for all inputs $v<2^{2n}+2^{n-1}-1$, it holds that $-2^n\le b-a<2^n$. This is done in two cases:

1. $a\in \{0,...,2^n\}$: This corresponds to $v\in \{2^{n-1},...,2^{2n}+2^{n-1}-2\}$.

   Two cases:

   1. $a\le b⟹0\le b-a<2^n-1⟹-2^n\le b-a<2^n⟹\delta =0$.
   2. $a>b⟹-2^n\le b-a<2^{n-1}⟹-2^n\le b-a<2^n⟹\delta =0$.

2. $a=-1$: This corresponds to $v\in \{-2^{n-1}+1,...,2^{n-1}-1\}$.

   $a=-1⟹1\le b-a<2^n⟹-2^n\le b-a<2^n⟹\delta =0$.

   Note: While this case algebraically includes negative values for $v$, the bounds $v\ge 0$ are implied by $v\in \mathbb{N}$. So this cannot be taken as proof that the approximation works for negative inputs.

Taken together, this proves $0\le v<2^{2n}+2^{n-1}-1⟹\delta =0$, which means that the approximation is exactly equal to rounded division by $2^n-1$ for those input.

Further, proving that the approximation is not exact for $v=2^{2n}+2^{n-1}-1$ is easy. This number corresponds to $a=2^n+1,b=0$, which results in a non-zero error term:

Therefore, $v=2^{2n}+2^{n-1}-1$ is the smallest (non-negative) input for which the approximation is not exact.

Showing failure points

After proving that the approximation fails, I thought it might also be interesting to see it fail for a few values of $n$. So here are the smallest inputs where the approximation starts to differ from the actual result for $n$ from 1 to 15:

These numbers were found experimentally using brute-force search.

Extending the input range

Depending on the use case, an input range of $v<2^{2n}+2^{n-1}-1$ might be too small. However, the trick can be extended to support larger inputs.

The main insight here is this equality:

Derived from:

This makes it possible to recursively rewrite the division by $2^n-1$ into a series of divisions by $2^n$.

Let $R_i$ be the approximation of rounded division by $2^n-1$ after $i$ iterations, defined as follows:

The trick from above then corresponds to $R_2$.

Let's see the smallest inputs $v$ the approximations $R_i$ start to fail for:

These numbers were found experimentally using a brute-force search. I aborted the search when it took too long, so some entries are missing.

The pattern is very clear. It seems that $R_i$ is exact for all inputs $v<2^{in}+2^{n-1}-1$. So the trick can be extended to support arbitrarily large input ranges by increasing the number of iterations.

In code, this is implemented as follows:

fn div_round_2pn_m1_iters(v: u32, n: u32, iters: u8) -> u32 {
    let round = 1 << (n - 1);
    let w = v + round;
    let mut r = w >> n; // R_1
    for _ in 1..iters {
        r = (r + w) >> n; // R_{i+1}
    }
    r
}

Proving correctness for extended ranges

Let $v,n,a,b$ be defined as before: $v=a(2^n-1)+b+2^{n-1}$. Further, let $i\in {\mathbb{N}}_1$ be the i-th iteration of the approximation $R_i$ as defined above and ${\delta }_i:=R_i-\mathrm{round}(v/(2^n-1))$. As before, $R_i$ is exact for inputs $v$ iff ${\delta }_i=0$.

I will start by simplifying the error term ${\delta }_i$:

With the error term in a nicer form, it's now easy to show by induction that $a=-1⟹{\delta }_i=0$:

• Base case ($i=1$): $a=-1⟹{\delta }_1=\lfloor (b-a)/2^n\rfloor =\lfloor (b+1)/2^n\rfloor =0$.
• Induction step: ${\delta }_{i+1}=\lfloor ({\delta }_i+b+1)/2^n\rfloor =\lfloor (b+1)/2^n\rfloor =0$.

Since $a=-1$ corresponds to $v\in \{-2^{n-1}+1,...,2^{n-1}-1\}$, this shows that the approximations (for any number of iterations) are exact for this range. (Again, the bounds $v\ge 0$ are implied by $v\in \mathbb{N}$.)

To make things simpler going forward, I derived an explicit formula for ${\delta }_i$ by unrolling the recursion:

Using the explicit formula for ${\delta }_i$, it's easy to show that $v=2^{in}+2^{n-1}-1$ results in ${\delta }_i=-1$, making the approximation non-exact for that input. $v=2^{in}+2^{n-1}-1$ corresponds to $a={\sum }_{l=0}^{i-1}{2}^{ln},b=0$. Plugging this into the explicit formula for ${\delta }_i$ gives:

Before I can prove the rest, I need split $a$ similar to how I split $v$ into $a,b$. Let $a=j_i{2}^{in}+k_i$ for $j_i\in \mathbb{Z},j_i:=\lfloor a/2^{in}\rfloor$ and $k_i\in \mathbb{N},k_i<2^{in},k_i:=a\mathrm{mod}{2}^{in}$ at iteration $i$.

Plugging this into the explicit formula for ${\delta }_i$ gives:

This way of representing ${\delta }_i$ reveals that ${\delta }_i$ is just $-j_i$ plus some small correction term ${\gamma }_i$ that depends on $k_i$ and $b$ but not $j_i$.

It's rather easy to show that ${\gamma }_i\in \{-1,0\}$ is always the case, which implies that ${\delta }_i\in \{-j_i-1,-j_i\}$.

However, that's not all that can be shown about ${\gamma }_i$. There is more structure to it. This structure becomes obvious by looking at a few concrete values. So are the possible values of ${\gamma }_2$ (with $n=2$) depending on $k_2$:

The table shows that ${\gamma }_2$ starts being always zero for $k_2<5$ and then switches to being either $-1$ or $0$ for $k_2\ge 5$.

This pattern persists across all values for $i,n$, but the switch occurs at different values $k_i$. I will call the value $k_i$ where the switch occurs the critical value, denoted as $c_i$. Here are a few values of $c_i$ I determined experimentally (missing values were too slow to compute):

The pattern is $c_i={\sum }_{l=0}^{i-1}{2}^{ln}=(2^{in}-1)/(2^n-1)$.

Showing that $k_i<c_i⟹{\gamma }_i=0$ is easy.

With this proven, 2 facts about ${\delta }_i$ are now known:

1. ${\delta }_i\in \{-j_i-1,-j_i\}$.
2. $k_i<c_i⟹{\delta }_i=-j_i$.

Since $v=2^{in}+2^{n-1}-1$ corresponds to $a={\sum }_{l=0}^{i-1}{2}^{ln}=c_i$ and $b=0$, it follows from (2) that $a\in \{0,...,c_i-1\}⟹{\delta }_i=0$. This range of $a$ corresponds to inputs $v\in \{2^{n-1},...,2^{in}+2^{n-1}-2\}$. Together with the earlier result for $a=-1$, this proves that the approximation $R_i$ is exact for all inputs $v<2^{in}+2^{n-1}-1$. $\square$

Other rounding modes

As it turns out, other rounding modes can also be implemented by making a small modification to the approximation. Similar to how $R_i$ was defined, $F_i$ and $C_i$ can be defined for floor and ceiling division, respectively:

where:

$F_i$ is exact for all inputs $v<2^{in}+2^n-2$, and $C_i$ is exact for all inputs $v<2^{in}$.

Proving correctness for floor division

The proof for floor division is very similar to the proof for rounded division. The main difference is that $v$ is split differently. Let $v=a_F(2^n-1)+b_F+2^n-1$ with $a_F,b_F$ defined similarly to before. Let ${\delta }_i^F:=F_i-\lfloor v/(2^n-1)\rfloor$ be the error term for floor division after $i$ iterations. Simplifying the error term gives:

Note that ${\delta }_i^F$ has the same recursive structure as ${\delta }_i$ (rounded division). Therefore, the rest of the proof follows the same steps as before, leading to the conclusion that $F_i$ is exact for all inputs $v<2^{in}+2^n-2$. Writing out this proof will be left as an exercise to the reader. $\square$

Proving correctness for ceiling division

Same game again. Let $v=a_C(2^n-1)+b_C+1$ with $a_C,b_C$ defined similarly to before. Let ${\delta }_i^C:=C_i-\lceil v/(2^n-1)\rceil$ be the error term for ceiling division after $i$ iterations. Simplifying the error term gives: