整數定數除法的代換 (constant integer division)

本文內容採用創用 CC 姓名標示-非商業性-相同方式分享 3.0 台灣 授權條款授權.

---

不論是 CPU 或是電路設計中，除法的速度都相當慢，是個無論如何都要避免的存在。一般而言，對於可以多次使用的除法，可以用計算倒數之後拿來乘（註：請見文尾）。而本文想要探討的是「整數的常數除法 (constant integer division)"」，像是「一個 32-bit 的 unsigned integer 除以 10」這種問題。

不妨來看一個實際例子，在 gcc 裡面有一個 32-bit 的 unsigned integer x，那麼 x/10 會被轉換成 (x*3435973837)>>35，可以透過這個指令觀察：

1echo "unsigned int f(unsigned int x){return x/10;}" | g++ -S -xc - -o -

在現代編譯器中，用乘法替換掉整數除法是非常成熟技巧，使用者自然也不須要知道其原理，但是知道背後的原理也是很有趣的，探討此一原理就是本篇文章的目的，而內容則是參考了 libgmp 大數運算 library 的文件。

實做的想法很簡單，既然答案是整數的話只要足夠近似就行了，然而「足夠」的定義是什麼呢？想當然 (x*3)>>5 或是 (x*26)>>8 都提供了類似 x/10 的功能，但是他們的結果都不完全正確：

這當然是不符合標準的，因此需要有更嚴謹的數學方法，接下來的文章中說明在除數、被除數都是 unsigned integer 的情形下如何找出最小的合適數字。

問題定義

首先我們都知道我們的被除數的範圍，所以我們的問題定義變成找出最小的 $(l, m)$ 使得（若存在的話）：

$$ \lfloor \frac{nm}{2^{N+l}} \rfloor = \lfloor \frac{n}{d} \rfloor, \quad \forall 0 \le n < 2^N$$

為了方便討論，我們多定義 $(q, r)$

$$ \left\lbrace \begin{aligned} q &= \lfloor \frac{n}{d} \rfloor \\ r &= n - qd \quad (0 \le r < d) \end{aligned} \right. $$

我們的問題定義可以改寫如下

$$ 0 \le \frac{nm}{2^{N+l}} - q < 1, \quad \forall 0 \le n < 2^N$$

仔細推導

先觀察左側的不等式，我們可以（小心地）這樣變化：

$$ \begin{aligned} & \quad 0 \le \frac{nm}{2^{N+l}} - q \\ \longrightarrow & \quad 0 \le \frac{nmd}{2^{N+l}} - dq \\ \longrightarrow & \quad 0 \le \frac{nmd}{2^{N+l}} - n + r \\ \longrightarrow & \quad 0 \le n\left(\frac{md}{2^{N+l}}-1\right) + r \end{aligned} $$

因為 $(n, r)$ 恆正，所以 $n$ 後面乘的那串 $\ge 0$ 就是等式成立的充分條件：

$$ \frac{md}{2^{N+l}}-1 \ge 0 \quad \longleftrightarrow \quad md \ge 2^{N+l} $$

右側的不等式類似：

$$ \begin{aligned} & \quad \frac{nm}{2^{N+l}} - q < 1 \\ \longrightarrow & \quad \frac{nm}{2^{N+l}} - qd - r < d-r \\ \longrightarrow & \quad n(\frac{md}{2^{N+l}} - 1) < d-r \end{aligned} $$

（可以發現證明中最重要的步驟就是把所有不等是的 $q$ 設法換成 $n$。）我們知道 $d-r \ge 1$，所以下面這個是充分條件：

$$ n(\frac{md}{2^{N+l}}-1) < 1$$

又，$n < 2^N$ 所以繼續放寬，下面這個也是充分條件：

$$ 2^N(\frac{md}{2^{N+l}}-1) \le 1 \quad \longleftrightarrow \quad md \le 2^{N+l}+2^N$$

結合兩個式子，我們需要的條件就是

$$ 2^{N+l} \le md \le 2^{N+l}+2^l $$

既然條件出來了，就可以透過暴力搜尋找出最小的 $(l, m)$，重新比對 gcc 把 x/10 變成 (x*3435973837)>>35 的步驟——在 $N = 32, d = 10$ 的前提下：

$l = 1 \quad \rightarrow \quad (2^{33}, 2^{34} + 2^1) = (8589934592, 8589934594)$，中間沒有 $d$ 的倍數。

$l = 2 \quad \rightarrow \quad (2^{34}, 2^{33} + 2^2) = (17179869184, 17179869188)$，中間也沒有 $d$ 的倍數。

$l = 3 \quad \rightarrow \quad (2^{35}, 2^{35} + 2^3) = (34359738368, 34359738376)$，中間終於有個 $d$ 的倍數了，因此 $(l, m) = (3, 3435973837)$。

文尾

還記得文章開頭有一句「對於可以多次使用的除法，可以用計算倒數之後拿來乘（註：請見文尾）」嗎，我想讀者應該是忘記了吧，因為我也差點忘記了。這個技巧最常用的地方是計算向量的 normalization。

$$ \frac{\mathbf{v}}{|\mathbf{v}|} = \mathbf{v} \times \frac{1}{|\mathbf{v}|}$$

本來左邊要算 $x$ 次除法，右邊只要一次除法跟 $x$ 次乘法。