注：本文中的算法來自于 Division by Invariant Integers using Multiplication ^[1]。

---

眾所周知，編譯器可以把變量除以常量優化為乘法和移位。 例如：

Uint32 f(Uint32 a) { return a / 3; }

會生成下面這樣的匯編（x86_64）：

f:mov eax, edimov edx, 2863311531imul rax, rdxshr rax, 33ret

你一定很好奇，這個 2863311531 是怎么得到的？在這之前，我們需要先準備一些預備知識。

首先我們需要知道，對于常見的指令集，都提供了N位整數乘以N位整數得到2N位結果的指令。以32位為例，對于 x86 指令集來說：

• mul r/m32：計算 eax 和 r/m32 進行無符號乘法的結果，將低32位存入 eax，高32位存入 edx。
• imul r/m32：計算 eax 和 r/m32 進行有符號乘法的結果，將低32位存入 eax，高32位存入 edx。

對于 ARM 指令集來說：

• umull RdLo, RdHi, Rm, Rs：計算 Rm 和 Rs 進行無符號乘法的結果，將低32位存入 RdLo，高32位存入 RdHi。
• smmul Rd, Rm, Rs：計算 Rm 和 Rs 進行有符號乘法的結果，將高32位存入 Rd。

我們可以用下面兩個函數代表計算乘法高位的結果：

Uint32 muluh(Uint32 a, Uint32 b) { return (Uint64(a) * b) >> 32; } Int32 mulsh(Int32 a, Int32 b) { return (Int64(a) * b) >> 32; }

這些指令為我們的優化提供了可能性。 為了充分利用這些指令，我們的目標是將 n / d 優化為 muluh(n, m) >> l，寫成數學表達式就是

。 那么，我們應該如何通過 計算出 和 呢？

定理1 設

，， 是非負整數，，且滿足 ，則對于 的所有整數 ，有 。

證明 設

，則由題設可得 。 對于 ，我們可以將 寫成 的形式，其中 ， 。 也就是說，我們需要證明 。

因為

， ， ，所以

所以

，得證。

換而言之，

的取值范圍是 。 我們可以先令 ，這樣可以保證取值范圍不為空，然后不斷減小 ，直到取值范圍內只有一個整數。 這樣可以使得 盡量小，并且在少數情況下可以使得 達到0，從而省略最后的移位。 （例如對于32位無符號乘法，除以641時就能省略最后的移位。）

代碼如下（注：本文中均以32位乘除法為例）：

constexpr int N = 32;inline int clz(Uint32 x) { return __builtin_clz(x); }struct Multiplier {Uint64 m;int l; };Multiplier chooseMultiplier(Uint32 d) {assert(d != 0);// l = ceil(log2(d))int l = N - clz(d - 1);Uint64 low = (Uint64(1) << (N + l)) / d;Uint64 high = ((Uint64(1) << (N + l)) + (Uint64(1) << l)) / d;while((low >> 1) < (high >> 1) && l > 0)low >>= 1, high >>= 1, --l;return {high, l}; }

試驗一下：

chooseMultiplier(3, 32);

得到

{2863311531, 1}

也就是說，我們可以將 n / 3 優化為

return muluh(n, 2863311531) >> 1;

然而，這個算法是存在問題的。例如，當

時，我們會得到這樣的結果：{4908534053, 3}

發現了嗎？

，超過了 Uint32 能夠表示的范圍！ 究其原因，是因為原始的范圍中就只有一個整數，無法縮小，而 。

不過不用擔心，注意到

，也就是說， 最大不會達到 。 所以我們可以利用乘法分配律將 拆成 和 兩個部分。 其中 ，可以像之前一樣用 muluh 解決，而 的高位結果直接就是 自己了。 也就是說：Uint32 t = muluh(n, m - (Uint64(1) << 32); return (n + t) >> l;

但是這樣還是有一點問題，如果

很大的話， 還是可能會溢出。 所以我們需要再進行一次變形：Uint32 t = muluh(n, m - (Uint64(1) << 32); return (((n - t) >> 1) + t) >> (l - 1);

也就是

。 將 代入的話，結果就是：Uint32 t = muluh(n, 613566757); return (((n - t) >> 1) + t) >> 2;

另外還有一種情況就是當

為偶數時，我們可以將 拆成 。 比如當 的時候，我們可以將它拆成 和 。 你可能會好奇， 時的 不是超過了 Uint32 的范圍嗎？ 實際上，在前一步的 之后，我們的除法實際上只需要31位的有效精度，在這個精度下計算出的 是不會超過范圍的。 加入有效精度限制后的 chooseMultiplier 如下：Multiplier chooseMultiplier(Uint32 d, int p) {assert(d != 0);assert(p >= 1 && p <= N);// l = ceil(log2(d))int l = N - clz(d - 1);Uint64 low = (Uint64(1) << (N + l)) / d;Uint64 high = ((Uint64(1) << (N + l)) + (Uint64(1) << (N + l - p))) / d;while((low >> 1) < (high >> 1) && l > 0)low >>= 1, high >>= 1, --l;return {high, l}; }

（需要注意的是，當

的時候，low 和 high 的計算過程中是會產生溢出的，但這種情況下除法的結果只可能是0或1，可以直接用比較解決，所以這里不做考慮。)

調用

chooseMultiplier(7, 31);

的結果是

{2454267027, 2}

所以 n / 14 可以被優化為：

return muluh(n >> 1, 2454267027) >> 2;

當然，眾所周知，如果

可以寫作 的形式的話，那么除法就可以直接優化為一個移位了。

另外需要注意的是，如果按完整精度計算出的

沒有達到 的話，就沒有必要進行這個操作，否則反而可能會使得結果更差。 例如對于 來說，直接計算的結果 muluh(n, 2863311531) >> 2 優于先除以二的結果 muluh(n >> 1, 2863311531) >> 1（后者多了一次移位）。

結合上述幾種情況，完整代碼如下：

#include <cassert> #include <initializer_list> #include <iostream>using Uint32 = unsigned int; using Uint64 = unsigned long long; using Int32 = int; using Int64 = long long;inline int clz(Uint32 x) { return __builtin_clz(x); } inline int ctz(Uint32 x) { return __builtin_ctz(x); }Uint32 muluh(Uint32 a, Uint32 b) { return (Uint64(a) * b) >> 32; } Int32 mulsh(Int32 a, Int32 b) { return (Int64(a) * b) >> 32; }constexpr int N = 32;struct Multiplier {Uint64 m;int l; };Multiplier chooseMultiplier(Uint32 d, int p) {assert(d != 0);assert(p >= 1 && p <= N);// l = ceil(log2(d))int l = N - clz(d - 1);Uint64 low = (Uint64(1) << (N + l)) / d;Uint64 high = ((Uint64(1) << (N + l)) + (Uint64(1) << (N + l - p))) / d;while((low >> 1) < (high >> 1) && l > 0)low >>= 1, high >>= 1, --l;return {high, l}; }void generateUnsignedDivision(Uint32 d) {assert(d != 0);std::cout << "Uint32 div" << d << "(Uint32 n) {n";if(d >= (Uint32(1) << (N - 1))) {std::cout << " return n >= " << d << ";n";} else {int s = ctz(d);if(d == (Uint32(1) << s)) {std::cout << " return n";if(s > 0) std::cout << " >> " << s;std::cout << ";n";} else {Multiplier multiplier = chooseMultiplier(d, N);if(multiplier.m < (Uint64(1) << N)) s = 0;else multiplier = chooseMultiplier(d >> s, N - s);if(multiplier.m < (Uint64(1) << N)) {std::cout << " return muluh(n";if(s > 0) std::cout << " >> " << s;std::cout << ", " << multiplier.m << ")";if(multiplier.l > 0) std::cout << " >> " << multiplier.l;std::cout << ";n";} else {std::cout << " Uint32 t = muluh(n, " << (multiplier.m - (Uint64(1) << N)) << ");n";std::cout << " return (((n - t) >> 1) + t) >> " << (multiplier.l - 1) << ";n";}}}std::cout << "}n"; }int main() {for(Uint32 d : std::initializer_list<Uint32>{1, 3, 6, 7, 14, 31, 32, 641, 0x7FFFFFFF, 0x80000000})generateUnsignedDivision(d); }

下面展示了該代碼生成的各種情況的優化結果：

Uint32 div1(Uint32 n) {return n; } Uint32 div3(Uint32 n) {return muluh(n, 2863311531) >> 1; } Uint32 div6(Uint32 n) {return muluh(n, 2863311531) >> 2; } Uint32 div7(Uint32 n) {Uint32 t = muluh(n, 613566757);return (((n - t) >> 1) + t) >> 2; } Uint32 div14(Uint32 n) {return muluh(n >> 1, 2454267027) >> 2; } Uint32 div31(Uint32 n) {Uint32 t = muluh(n, 138547333);return (((n - t) >> 1) + t) >> 4; } Uint32 div32(Uint32 n) {return n >> 5; } Uint32 div641(Uint32 n) {return muluh(n, 6700417); } Uint32 div2147483647(Uint32 n) {Uint32 t = muluh(n, 3);return (((n - t) >> 1) + t) >> 30; } Uint32 div2147483648(Uint32 n) {return n >= 2147483648; }

---

這回講了無符號除法的優化，關于有符號除法，我們下次再說。

參考

^Torbj?rn Granlund and Peter L. Montgomery. 1994. Division by invariant integers using multiplication. SIGPLAN Not. 29, 6 (June 1994), 61–72.?https://dl.acm.org/doi/10.1145/773473.178249

總結

以上是生活随笔為你收集整理的无符号有符号乘法_【编译笔记】变量除以常量的优化（一）——无符号除法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。