Linux 核心原始程式碼的整數除法

--- tags: LINUX KERNEL, LKI --- # Linux 核心原始程式碼的整數除法 > 資料整理: [jserv](http://wiki.csie.ncku.edu.tw/User/jserv), [jouae](https://github.com/jouae) ## 簡介 Linux 核心原始程式碼蘊含多項巧思，本文討論其針對整數除法的實作手法。 ## 整數除法的向上取整在數學和電腦科學中，[取整函數](https://en.wikipedia.org/wiki/Rounding)的作用將實數映射到相近整數，常見手法有「下取整數」(floor)和「上取整數」(ceiling)。其中，上取整函數即為「取頂符號」，在數學記作 $\lceil x \rceil$，在電腦科學記作 `ceil(x)`，表示不小於 $x$ 的整數中最小的數值。例如：$\lceil 2.718281 \rceil = 3$。為了避免後續討論出現歧意，我們先回顧 C 語言在除法運算子 `/` 對整數型態的數值運算定義。C99 規格書中 6.5.5 Multiplicative operators 第 6 點指出 C 語言中除法 `/` 運算子，在兩整數相除後，其小數部分是直接捨去。這個做法也被稱為截斷至 $0$ (truncation toward zero)。例如：$7/(-3)=-2$ 與 $7/3=2$。 > When integers are divided, the result of the / operator is the algebraic quotient with any fractional part discarded.^87^) (87) This is often called "truncation toward zero" 為了避免後續討論混淆一般數學中的除法 (結果可能有餘數) 與該 C 語言除法運算子作用 (結果為整數)，在本文強調以下用語： * 除法：泛指數學中除法運算 * 整數除法：特指 C 語言中兩整數型態數值經由 `/` 運算子計算則二個整數 $n$ 和 $d$ 用文字表示數學除法或是除法運算子作用寫為： | 數學除法 | 除法運算子 | | -------- | -------- | | $n$ 除以 $d$ | $n$ 整數除以 $d$ | | $d$ 除 $n$ | $d$ 整數除 $n$ | 此處整數除並非整除，使用「整數除」一詞在此僅用於明確表示與數學除法的差異，並非用於解釋兩數其一為另一數的因數。 > [「除以」的意思](https://www.facebook.com/thinkinggarden/posts/pfbid027N51doVuhNksHGHHr7wekk632vYTP72QXEczG9DB9yKdDFdFvZHtjumpsGVqzpfFl) 給定二個整數 $n$ 和 $d$，如何將 $n$ 除以 $d$ 後進行向上取整呢？直觀的作法如下： ```c int div_and_ceil(int n, int d) { if (n % d == 0) return n / d; return n / d + 1; } ``` 但在 Linux 核心原始程式碼中，相同目的的程式碼會用以下形式的巨集呈現： (檔案: [include/linux/math.h](https://github.com/torvalds/linux/blob/master/include/linux/math.h) 和 [include/uapi/linux/const.h](https://github.com/torvalds/linux/blob/master/include/uapi/linux/const.h)) ```c #define DIV_ROUND_UP __KERNEL_DIV_ROUND_UP #define __KERNEL_DIV_ROUND_UP(n, d) (((n) + (d) - 1) / (d)) ``` 該巨集的使用前提是，$n$ 和 $d$ 都是正整數，否則會失敗。該手法結合除法和向上取整運算的數學原理與 C 語言整數除法的處理邏輯。接下來解析該操作：對一個數值進行上取整會得到什麼？當數字有非零的小數部分時，會得到下一個整數；否則，得到相同的數字。將 $n$ 除以 $d$ 會得到一個非零小數結果，當且僅當被除數 $n$ 不是除數 $d$ 的整數倍；換言之，餘數不為 `0`。關於說明上取整函式、餘數和商數 (quotient) 之間的關係，可見以下案例： ![image](https://hackmd.io/_uploads/H1_r6HWL0.png) 觀察上方案例得知，將除法結果做為上取整函數的輸入值，當餘數為非 $0$ 時，需要對商數加 $1$。為此，我們需要在被除數上加某個數值，使餘數為 $0$ 時，商數不變，而餘數為非 $0$ 時，商數增加 $1$，該數值就是 $d - 1$。以下觀察當我們將 $d - 1$ 加到被除數 $n$ 的過程。由於餘數 $r$ 永遠不會大於除數 $d$，亦即 $$ 0 \leq r < d $$ 加上 $d - 1$ 後，得到新的餘數 $r'$： $$ d-1 \leq r' < 2d-1 $$ 考慮 $r'$ 有以下二種狀況。 - [ ] 情況 1 當 $r$ 為 $0$ 時，加上 $d - 1$ 後，新的餘數 $r'$ 變成 $d - 1$。由於這是整數除法，$\dfrac{d-1}{d}$ 為 $0$，因此商數不變。 - [ ] 情況 2 當 $r$ 為非 $0$ 時，由於 $r$ 至少為 $1$，加入 $d - 1$ 使 $r'$ 至少為 $d$，最多為 $2d - 1$，這小於 $2d$。所以，加入 $d - 1$ 剛好在商數上加 $1$，因為 $\dfrac{r'}{d} = 1$。從[歐幾里德除法](https://en.wikipedia.org/wiki/Euclidean_division)思考，假設 $n$ 為被除數，$d$ 為除數同時不為 $0$，且 $n,d$ 皆為整數，則歐幾里德除法告訴我們，存在唯一的二整數，其一為商數 $q$ ，其一為整數餘數 $r$： $$ n=dq+r, \text{ and } 0 \leq r < \vert d \vert $$ 由於該巨集的使用限制 $d$ 固定為正數，所以不等式寫為 $0 \leq r < d$ 現加入一整數 $d-1$ 至等式兩側則 $$ n+(d-1) = dq+r +(d-1) = dq + (r+d-1) $$ 對等式兩側除以 $d$ $$ \dfrac{n+(d-1)}{d} = q + \dfrac{r+d-1}{d} $$ :warning: 第二項除式為何不整理成 $1+\dfrac{r-1}{d}$ 呢？此處的「除法」其實是 C 語言中的「整數除法」，倘若簡寫成 $1+\dfrac{r-1}{d}$，C 語言「整數除法」會有非預期的結果，例如：當 $r=0$ 時，該簡式可寫為 $1+\dfrac{-1}{d}$，則該簡式中第二項經由整數除法得到該簡式結果為 $1$。使用原式計算的話，$\dfrac{r+d-1}{d} = \dfrac{d-1}{d}$ 經由整數除法得到該原式結果為 $0$。如上段討論情況 1 與情況 2 ，對等式右側第二項除式 $\frac{r+d-1}{d}$ 討論以下二種情況： - [ ] 情況 1 若 $r$ 為 $0$，則第二項除式為 $\dfrac{d-1}{d}$，在 C 語言整數除法中，該整數除式 $\dfrac{d-1}{d}$ 會被截斷至 $0$。 - [ ] 情況 2 若 $r$ 不為 $0$，因為 $r$ 為整數，則 $r$ 至少為 $1$，藉由該整數餘數不等式 $1 \leq r < d$ 第二項除式滿足不等式 $d \leq r + d-1 < 2d - 1$，於是在 C 語言整數除法中，該整數除式 $\dfrac{r+d-1}{d}$ 會被截斷至 $1$。對應的圖解： ![image](https://hackmd.io/_uploads/BJi91UbIA.png) 因此，$d - 1$ 是恰到好處的數值，它能在餘數為 $0$ 時，使商數保持不變，並在餘數為非 $0$ 時，將商數增加 $1$。這就是 `DIV_ROUND_UP` 的原理，Linux 核心開發者改寫二個分支為不需要分支的精簡實作，展現數學之美。至於整數除法的四捨五入，可參見 `DIV_ROUND_CLOSEST` 巨集 (檔案: [include/linux/math.h](https://github.com/torvalds/linux/blob/master/include/linux/math.h)): ```c /* Divide positive or negative dividend by positive or negative divisor * and round to closest integer. Result is undefined for negative * divisors if the dividend variable type is unsigned and for negative * dividends if the divisor variable type is unsigned. */ #define DIV_ROUND_CLOSEST(x, divisor)({ \ typeof(x) __x = x; \ typeof(divisor) __d = divisor; \ (((typeof(x))-1) > 0 || \ ((typeof(divisor))-1) > 0 || \ (((__x) > 0) == ((__d) > 0))) ? \ (((__x) + ((__d) / 2)) / (__d)) : \ (((__x) - ((__d) / 2)) / (__d)); \ }) ``` `DIV_ROUND_CLOSEST` 巨集執行除法並將結果四捨五入到最接近的整數，可處理正負數的被除數和除數。原理如下： 1. 檢查被除數或除數是否為無號整數。若是，或若被除數和除數都為正或負，則將除數的一半加到被除數，再除以除數 2. 若被除數和除數的正負號不同，則從被除數中減去除數的一半，再除以除數該巨集使用 $\dfrac{dividend + \frac{divisor}{2}}{divisor}$ 來對正數結果進行四捨五入，並避免在正負號不同和無號整數類型時出現錯誤結果。 ## `do_div` 巨集 `do_div` 是 Linux 核心原始程式碼裡頭的巨集，進行整數除法操作，將商數儲存在給定的變數中，並將餘數儲存在另一個變數。巨集定義如下 (檔案: [include/asm-generic/div64.h](https://github.com/torvalds/linux/blob/master/include/asm-generic/div64.h)) ```c # define do_div(n,base) ({ \ uint32_t __base = (base); \ uint32_t __rem; \ __rem = ((uint64_t)(n)) % __base; \ (n) = ((uint64_t)(n)) / __base; \ __rem; \ }) ``` 參數說明： - `n`：要進行除法的 64 位元整數。 - `base`：除數，一個 32 位元整數。巨集的返回值是 $n$ 除以 $base$ 的餘數。用法： ```c u64 num = 1234567890123456ULL; u32 base = 1000; u32 remainder; remainder = do_div(num, base); ``` `num` 被 1000 除，商數則儲存在 `num` 中，餘數儲存在 `remainder` 中。使用 `do_div` 巨集的好處是，可一次得到商數和餘數，且運算過程不依賴暫存變數，更重要的是，`do_div` 巨集可消除非預期的連結錯誤：針對特定的 Arm 處理器，某些 GCC Toolchain 會將 (當被除數是) 64 位元除法轉換為 `__aeabi_uldivmod` 呼叫 (這是 [libgcc](https://gcc.gnu.org/onlinedocs/gccint/Libgcc.html) 的一部分)，即使除數是 32 位元，這仍然會執行不必要的 64 位元的除法操作。由於效率原因，Linux 核心中不實作這類 64 位元除法，導致連結時出現未定義符號的引用錯誤。因此，對 64 位元除法應使用 `do_div` 巨集。 > [[PATCH] sched: Explicit division calls on 64-bit integers](https://lists.infradead.org/pipermail/linux-arm-kernel/2012-November/133312.html) 注意：上述巨集定義是針對 64 位元的處理器架構，而針對 32 位元的處理器架構，Linux 核心另行提供其他的實作方式，以避免非預期的效能議題，可參見 [[PATCH 2/5] do_div(): generic optimization for constant divisor on 32-bit machines](https://lore.kernel.org/lkml/1446503610-6942-3-git-send-email-nicolas.pitre@linaro.org/)。 `do_div()` 巨集針對 Linux 核心原始程式碼，後者幾乎不使用 FPU 指令，但正如 Daniel Lemire 教授在 [Fast exact integer divisions using floating-point operations](https://lemire.me/blog/2017/11/16/fast-exact-integer-divisions-using-floating-point-operations/) 一文指出，整數除法改用 FPU 運算，可帶來更好的效能表現，也就是說，`do_div()` 巨集不見得適用於使用者層級的程式碼。 > [Re: do_div considered harmful](https://lkml.iu.edu/hypermail/linux/kernel/0308.0/0617.html) ## 針對除數已知狀況的加速除法運算的成本往往高於加法和乘法，特別在某些硬體平台缺乏整數除法指令 (如 Arm Cortex-A8) 的狀況，就不得不仰賴高成本的除法的模擬運算，而即使 Intel 公司出品的 Core 系列處裡器 [Sandy Bridge 微架構](https://en.wikipedia.org/wiki/Sandy_Bridge)中，針對 64 位元整數的除法運算，會佔用 40 到 103 個處理器週期，運算成本仍屬沈重。 > 來源: [Agner Fog’s instruction tables](https://www.agner.org/optimize/instruction_tables.pdf)，第 180 頁既然除法運算的成本居高不下，於是工程人員就著手改進整數除法的策略。其一是將除法改為乘法運算，例如 $\frac{n}{d}$ 在分子和分母分別乘上 $2^N$ 後，得到等價的 $n\times\frac{\frac{2^N}{d}}{2^N}$，其中對 2 的 N 次方 (power of 2) 進行除法，在無號整數的操作就等同於右移 (shift) N 個位元，又當 $d$ (除數) 的數值是已知的狀況下，$\frac{2^N}{d}$可預先計算，也就是說，$\frac{n}{d}$ 可轉換為乘法和位移操作，進而減少運算成本。不過，顯然 $\frac{2^N}{d}$ 無法總是用整數來表達 (僅在 $d$ 是 2 的冪時才有效)，因此，我們可先估算，並將差值補償回去。重新檢視 $\frac{n}{d} = n \times \frac{\frac{2^N}{d}}{2^N}$，當我們想將 $n$ 除以 $4$ 時，就相當於乘以 $0.25$，所以我們可將 $\frac{5}{4}$ 改為 $5 \times 0.25$，我們得到整數的部分 (即 $1$)，和小數部分 (即 $0.25$)，後者乘以 $4$ 即是 $1$，也就是餘數。下方程式碼展示這概念: ```c const uint32_t D = 3; #define M ((uint64_t)(UINT64_C(0xFFFFFFFFFFFFFFFF) / (D) + 1)) /* compute (n mod d) given precomputed M */ uint32_t quickmod(uint32_t n) { uint64_t quotient = ((__uint128_t) M * n) >> 64; return n - quotient * D; } ``` 其中 `__uint128_t` 是 GCC extension，用以表示 128 位元的無號數值。 > [GCC: C Extensions: 128-bit Integers](https://gcc.gnu.org/onlinedocs/gcc/_005f_005fint128.html) $n$ 模除 $3$ 的數學推導如下： $$ \begin{align} n \% 3 &= n - \left \lfloor {\frac{n}{3}} \right\rfloor \times 3 \\ &= n - \left \lfloor n \times \frac{ \frac{2^{64}}{3} }{2^{64}} \right\rfloor \times 3\\ &= n - \left \lfloor \left( n \times \frac{2^{64}}{3} \right) \times \frac{1}{2^{64}} \right \rfloor \times 3\\ &= n - ((n \times M) >> 64) \times 3 \\ \end{align} $$ 其中 $\frac{2^{64}}{3}$ 對應到 `M` 無條件捨去的除法運算再補 $1$ 對應到上取整數，故 $M = \left\lceil \frac{2^{64}}{3} \right \rceil$ 由 Meta 公司維護的 [jemalloc](https://github.com/jemalloc/jemalloc) 是個高效的記憶體配置器 (memory allocator，即 malloc 和 free 函式對應的實作)，特別在多核多執行緒的環境有優異表現，在其原始程式碼 [include/jemalloc/internal/div.h](https://github.com/jemalloc/jemalloc/blob/dev/include/jemalloc/internal/div.h) 和 [src/div.c](https://github.com/jemalloc/jemalloc/blob/dev/src/div.c) 也運用類似的技巧：已知 $\frac{n}{d} = n \times \frac{\frac{2^k}{d}}{2^k}$ 但 $k$ 要取多少，方可滿足我們要的準確度呢? [src/div.c](https://github.com/jemalloc/jemalloc/blob/dev/src/div.c) 提供證明。假設 $n, d, \frac{n}{d}$ 為整數 $$ \begin{align*} \frac{n}{d} &= p, \ p \in Z \\ &= \left \lfloor \left\lceil \frac{2^k}{d}\right\rceil \times \frac{n}{2^k}\right \rfloor , \ \forall k \in N\\ &= \left \lfloor \frac{2^k + r}{d} \times \frac{n}{2^k}\right \rfloor (\because r = - 2^k \mod d)\\ &= \left \lfloor \frac{n}{d} + \frac{r}{d} \times \frac{n}{2^k} \right \rfloor \\ &= \frac{n}{d} + \left \lfloor \frac{r}{d} \times \frac{n}{2^k} \right \rfloor (\because \frac{n}{d} = p \in Z) \\ \end{align*} $$ 證明 $\left\lceil \frac{2^k}{d}\right\rceil = \frac{2^k + r}{d}$ $$ \left\lceil \frac{2^k}{d}\right\rceil = \frac{2^k + d-(2^k \mod d)}{d} = \frac{2^k + r}{d} =\frac{2^k + (-2^k \mod d)}{d} $$ 證明 $-2^k \mod d= d-(2^k \mod d)$ $$ \begin{align*} \Rightarrow & 2^k \mod d = r \\ 2^k &= pd + r \\ -2^k &= (-p)d + -r \\ -2^k &= (-p-1)d + d-r \\ -2^k &= p'd + (d-r) \\ \Rightarrow & (-2^k \mod d ) = d-r\\ \end{align*} $$ $\frac{n}{d} = \frac{n}{d} + \left \lfloor \frac{r}{d} \times \frac{n}{2^k} \right \rfloor$ 成立的條件為 $\frac{r}{d} \times \frac{n}{2^k} < 1$ * 因為 $r = 2^k \mod d$，所以 $\frac{r}{d} <１$ * 如果 $\frac{n}{2^k} < 1$ 則 $\frac{r}{d} \times \frac{n}{2^k} < 1$ 恆成立 * 假設 n 的範圍是 unsigned int 則 k 取 32 可以保證 $\frac{n}{2^k} <1$ 以下是改寫 jemalloc 快速除法的精簡實作: ```c #include <limits.h> #include <stdbool.h> #include <stddef.h> #include <stdint.h> #include <stdio.h> #include <stdlib.h> /* find last (most-significant) bit set */ #define fls32(x) ((x) ? sizeof(int) * 8 - __builtin_clz(x) : 0) void fastdiv_prepare(uint32_t div, uint32_t *m, uint8_t *s1, uint8_t *s2) { const int l = fls32(div - 1); const uint64_t nextpow2 = UINT64_C(1) << l; uint32_t magic = ((nextpow2 - div) << 32) / div; *m = magic + 1; *s1 = (l > 1) ? 1 : l; *s2 = (l == 0) ? 0 : l - 1; if (div == 1) return; if (magic == 0) { /* div == nextpow2 */ *s1 -= 1, *s2 += 1; } else { magic = (magic + UINT64_C(0x100000000)) / 2; uint32_t rem = (nextpow2 << 31) - magic * div; assert(rem > 0 && rem < div); if (rem > div - (nextpow2 >> 1)) { *m = magic + 1; *s1 -= 1; } } } static uint32_t fastdiv(uint32_t val, uint64_t mul, uint8_t s1, uint8_t s2) { const uint32_t hi = (val * mul) >> 32; return (hi + ((val - hi) >> s1)) >> s2; } ``` 使用方式: ```c uint32_t div, mul; uint8_t s1, s2; div = 137; fastdiv_prepare(div, &mul, &s1, &s2); assert(fastdiv(1234, mul, s1, s2) == 1234 / div); ``` 比較 jemalloc 的快速除法與通用除法運算 * 編譯器最佳化程度由 `-O0` 到 `-O3` * 測量 10000 次畫出時間分佈圖測試結果： - [ ] 編譯器選項 `-O0` ![](https://hackmd.io/_uploads/r18-kqte2.png) - [ ] 編譯器選項 `-O1` ![](https://hackmd.io/_uploads/HyzGk5Flh.png) - [ ] 編譯器選項 `-O2` ![](https://hackmd.io/_uploads/HyRM1qYgn.png) - [ ] 編譯器選項 `-O3` ![](https://hackmd.io/_uploads/HyzV15Feh.png) 可見 jemalloc 的快速除法的確有效益。 ## 特別的除數儘管電腦採用二進位，但由於需要讓執行結果被人類所識別，因此不乏會有 mod 10 和 div 10 等操作，我們可考慮以下在單一函式取得商數和餘數、不依賴除法和乘法指令的實作： ```c #include <stdint.h> void divmod_10(uint32_t in, uint32_t *div, uint32_t *mod) { uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */ uint32_t q = (x >> 4) + x; x = q; q = (q >> 8) + x; q = (q >> 8) + x; q = (q >> 8) + x; q = (q >> 8) + x; *div = (q >> 3); *mod = in - ((q & ~0x7) + (*div << 1)); } ``` 以下先假設 $in \leq \cfrac{128 \times64}{67}\approx122$，隨後我們會看到為何有這個數值，並推廣到 $in \ge 122$。首先 `uint32_t x = (in | 1) - (in >> 2);` 等價於 $$ x= \begin{cases} \cfrac{3(in + 1)}{4} , & \text{if x is even} \\ \cfrac{3in}{4} , & \text{if x is odd} \end{cases} $$ 經過 `uint32_t q = (x >> 4) + x;` 後，等價於 $$ q= \begin{cases} \cfrac{67(in + 1)}{64} , & \text{if x is even} \\ \cfrac{67in}{64} , & \text{if x is odd} \end{cases} $$ 因為已假設 $in \leq 122$，所有的 `q = (q >> 8) + x;` 都可忽略；這是因右移 8 位元，也就是除以 128 後會為 $0$，僅剩 `+ x` 有值，故可以直接略過該四行敘述，令 `q` 等於上式。為了簡化說明，假設 `q` 為偶數；那麼右移 3 位元，也就是除以 8 後 $$ \cfrac{q}{8}=(\cfrac{3}{4}in+\cfrac{19}{64}in)\cdot\cfrac{1}{8}=\cfrac{3}{32}in+\cfrac{19}{512}in \approx div $$ 由於 $in < 122$，後項可忽略，得前項$\cfrac{3}{32}\approx\cfrac{1}{10}$，即是 `div`；而由於 $in=10div+mod$ ，可得 $$ mod = in - 10div \approx in- \cfrac{30}{32}in=in-(\text{(q & ~0x7)}+\cfrac{6}{32}in)=in-(\cfrac{24}{32}in+\cfrac{6}{32}in) $$ `(q & ~0x7)`可以視作右移 3 位後再左移 3 位，效果如同 $in-\cfrac{8}{32}in=\cfrac{24}{32}in$。當 $128+122 > in \geq122$ 時，`q = (q >> 8) + x;` 會被觸發，使得其值域再次回到 $in < 122$，故上述依然成立。 Linux 核心的二進位到十進位轉換程式碼，使用一些技巧，包括針對受限範圍使用較小的常數，以及一個 200 位元組的 mod-100-to-2-digits 查表方式。以下是相關程式碼： (部分刪減，檔案: [lib/vsprintf.c](https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/linux.git/tree/lib/vsprintf.c)) ```c /* Decimal conversion is by far the most typical, and is used for * /proc and /sys data. This directly impacts e.g. top performance * with many processes running. We optimize it for speed by emitting * two characters at a time, using a 200 byte lookup table. This * roughly halves the number of multiplications compared to computing * the digits one at a time. Implementation strongly inspired by the * previous version, which in turn used ideas described at * <http://www.cs.uiowa.edu/~jones/bcd/divide.html> (with permission * from the author, Douglas W. Jones). * * It turns out there is precisely one 26 bit fixed-point * approximation a of 64/100 for which x/100 == (x * (u64)a) >> 32 * holds for all x in [0, 10^8-1], namely a = 0x28f5c29. The actual * range happens to be somewhat larger (x <= 1073741898), but that's * irrelevant for our purpose. * * For dividing a number in the range [10^4, 10^6-1] by 100, we still * need a 32x32->64 bit multiply, so we simply use the same constant. * * For dividing a number in the range [100, 10^4-1] by 100, there are * several options. The simplest is (x * 0x147b) >> 19, which is valid * for all x <= 43698. */ static const u16 decpair[100] = { #define _(x) (__force u16) cpu_to_le16(((x % 10) | ((x / 10) << 8)) + 0x3030) _( 0), _( 1), _( 2), _( 3), _( 4), _( 5), _( 6), _( 7), _( 8), _( 9), _(10), _(11), _(12), _(13), _(14), _(15), _(16), _(17), _(18), _(19), _(20), _(21), _(22), _(23), _(24), _(25), _(26), _(27), _(28), _(29), ... _(90), _(91), _(92), _(93), _(94), _(95), _(96), _(97), _(98), _(99), #undef _ }; /* This will print a single '0' even if r == 0, since we would * immediately jump to out_r where two 0s would be written but only * one of them accounted for in buf. This is needed by ip4_string * below. All other callers pass a non-zero value of r. */ static char *put_dec_trunc8(char *buf, unsigned r) { unsigned q; /* 1 <= r < 10^8 */ if (r < 100) goto out_r; /* 100 <= r < 10^8 */ q = (r * (u64)0x28f5c29) >> 32; *((u16 *)buf) = decpair[r - 100*q]; buf += 2; /* 1 <= q < 10^6 */ if (q < 100) goto out_q; /* 100 <= q < 10^6 */ r = (q * (u64)0x28f5c29) >> 32; *((u16 *)buf) = decpair[q - 100*r]; buf += 2; /* 1 <= r < 10^4 */ if (r < 100) goto out_r; /* 100 <= r < 10^4 */ q = (r * 0x147b) >> 19; *((u16 *)buf) = decpair[r - 100*q]; buf += 2; out_q: /* 1 <= q < 100 */ r = q; out_r: /* 1 <= r < 100 */ *((u16 *)buf) = decpair[r]; buf += r < 10 ? 1 : 2; return buf; } ``` > 以下內容部分取自 [Ideal divisors: when a division compiles down to just a multiplication](https://lemire.me/blog/2021/04/28/ideal-divisors-when-a-division-compiles-down-to-just-a-multiplication/) 既然除法指令是 CPU 最耗時的指令之一。因此，編譯器最佳化的手段就是將已知常數的整數除法，改寫為乘法和位移。然而，有時編譯器甚至不需要位移，這樣的除數在本文稱為理想除數，與[費馬數](https://en.wikipedia.org/wiki/Fermat_number)有關。對於 32 位元無號整數，考慮二個理想除數 $641$ 和 $6700417$。對於 64 位元無號整數，考慮二個不同的理想除數$274177$ 和 $67280421310721$，後者分別是 $2^{32} + 1$ 和 $2^{64} + 1$ 的因數，且都是質數。於是： $$ \lfloor n / 274177 \rfloor = (n \times 67280421310721) \gg 64 $$ 和 $$ \lfloor n / 67280421310721 \rfloor = (n \times 274177) \gg 64 $$ 在這些表達式中，乘法是完全乘法（得到 128 位元的結果）。看似有 64 位元的位移，但實際上編譯後，位移就消失。不是所有編譯器都能做到這點，不過 GCC 和 Clang 可以。以下是 GCC 在 x86-64 平台上編譯 `n / 274177` 和 `n / 67280421310721` 時產生的組合語言輸出： ```c movabs rdx, 67280421310721 mov rax, rdi mul rdx mov rax, rdx ret mov rax, rdi mov edx, 274177 mul rdx mov rax, rdx ret ``` 在 Arm64 平台上也有類似的結果: ```c mov x1, 53505 movk x1, 0xf19c, lsl 16 movk x1, 0x3d30, lsl 32 umulh x0, x0, x1 ret mov x1, 12033 movk x1, 0x4, lsl 16 umulh x0, x0, x1 ret ``` 那麼餘數呢？優秀的編譯器會先計算商數，再藉由乘法和減法來推導餘數，不過在許多情況下，計算餘數比計算商數的代價更高。某些情況下，你可更高效計算餘數。論文〈[Faster Remainder by Direct Computation](https://arxiv.org/abs/1902.01961)〉提到額外的技巧：直接計算餘數，只要二個乘法： $$ n \% 274177 = \left( \text{UINT64_C}(n \times 67280421310721) \times 274177 \right) \gg 64 $$ 和 $$ n \% 67280421310721 = \left( \text{UINT64_C}(n \times 274177) \times 67280421310721 \right) \gg 64 $$ 換句話說，以下二個 C 函式完全等效： ```c // computes n % 274177 uint64_t div1(uint64_t n) { return n % 274177; } // computes n % 274177 uint64_t div2(uint64_t n) { return (UINT64_C(n * 67280421310721) * ((__uint128_t) 274177)) >> 64; } ``` 儘管第二個函式看似較冗長且複雜，但它通常會編譯成更高效的機械碼，只涉及二個乘法，通常比編譯器產生的最佳化指令更高效。延伸閱讀: [Integer Division by Constants: Optimal Bounds](https://arxiv.org/abs/2012.12369) / [fast_division](https://github.com/lemire/fast_division)