# Greatest Common Divisor 特性和實作考量 :::success 探討最大公因數 (GCD) 特性，考慮微處理器架構實作帶來的效能衝擊。考慮到 binary GCD 及其最佳化。 ::: Greatest Common Divisor 中譯「最大公因數」，簡稱 GCD。 給定整數 $a, b$，它倆各自有**自己的因數**[^factor]，取**相同**的因數中**最大**的數，即最大公因數 $\gcd(a, b)$ ## GCD 性質 - $\gcd(a, b) = \gcd(b, a)$ - $\gcd(a, 0) = |a|$ - $\gcd(a, b) = \gcd(a, b \% a)$ > $\%$ 是取餘數運算，表示存在 $k$ 滿足 $0 \leq b\% a = b - k\cdot a < a$ ## Euclidean algorithm > 歐幾里得法，也就是「輾轉相除法」 :::info 給定整數 $a, b$ 求出 $\gcd(a, b)$ ::: 令 $a_0=a, b_0=b$，根據上述 GCD 性質： $$ \begin{split} \gcd(a_0, b_0) &= \gcd(b_0 \% a_0 = \underline{b_1}, a_0) \\ \gcd(\underline{b_1}, a_0) &= \gcd(a_0 \% b_1 = \underline{a_1}, b_1) \\ \gcd(\underline{a_1}, b_1) &= \gcd(b_1 \% a_1 = b_2, a_1) \\ &\vdots \\ \gcd(a_i, b_i) &= \gcd(b_i \% a_i = a_{i+1}, b_i) \\ &\vdots \\ \gcd(\mathbf{0}, b_n) &= |b_n| \end{split} $$ ![](https://i.imgur.com/HtQU22S.gif) [^euclidean] > 展示動畫 綜合上面過程，可得程式碼實作： ```cpp int gcd(int a, int b) { return a ? gcd(b % a, a) : b; } ``` ## Bezout's Theorem 對於所有整數 $a, b$，存在**整數** $x, y$ 使得 $ax + by = \gcd(a, b)$ ## Extended Euclidean algorithm :::info 給定整數 $a, b$ 求出整數 $x, y$ 滿足 $ax + by = \gcd(a, b)$ ::: 令 $g = \gcd(a, b)$，配合上述 Bezout's Theorem 可得： $$ 0\cdot x_n + g\cdot y_n = g $$ 此處 $x_n \in \mathbb{Z}, y_n = 1$ > 即 $x_n$ 為任意整數 而根據 GCD 性質 $\gcd(a_i, b_i) = \gcd(b_i\%a_i, a_i)$ 以及 $b_i \% a_i = b_i - \lfloor{b_i\over a_i}\rfloor \cdot a_i$ 得到 $$ \begin{split} g &= x_j\cdot (b_i \% a_i) &+ y_j &\cdot a_i \\ &= x_j\cdot (b_i - \lfloor{b_i\over a_i}\rfloor \cdot a_i) &+ y_j &\cdot a_i \\ &= x_j \cdot b_i &+ (y_j - \lfloor{b_i\over a_i}\rfloor \cdot x_j) &\cdot a_i \end{split} $$ 也就是說 $$ g = a_i \cdot x_{j-1} + b_i \cdot y_{j-1} \text{ for }\begin{cases} x_{j-1} = y_j - \lfloor{b_i\over a_i}\rfloor \cdot x_j \\ y_{j-1} = x_j \end{cases} $$ 綜合上述過程： ```cpp /* It uses pointers to return multiple values */ int extended_gcd(int a, int b, int *x, int *y) { if (a == 0) { *x = 0; *y = 1; return b; } int _x, _y; int gcd = extended_gcd(b % a, a, &_x, &_y); *x = _y - (b/a) * _x; *y = _x; return gcd; } ``` 延伸閱讀: * [Extended Euclidean algorithm](https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm) ## GCD 實作和微處理器指令 回顧歐幾里得法: 1. 大數 := 大數 - 小數 2. 若「最大」和「最小」兩數相同，那麼它就是最大公因數，否則，回到第 1 步 示範: GCD(48, 40) ``` 48 - 40, 40 -> 8, 40 8 != 40, 於是重複上述動作: 40 - 8, 8 -> 32, 8 8 != 32, 於是重複上述動作: 32 - 8, 8 -> 24, 8 8 != 24, 於是重複上述動作: 24 - 8, 8 -> 16, 8 8 != 16, 於是重複上述動作: 16 - 8, 8 -> 8, 8 8 == 8 成立，於是 8 為 GCD ``` 用 C 語言實現歐幾里得法: (v1) ```cpp int findGCD(int a, int b) { while (1) { if (a > b) a -= b; else if (a < b) b -= a; return a; } } ``` ### 針對 [Pentium 4](https://en.wikipedia.org/wiki/Pentium_4) 硬體的實驗 以 Intel C compiler 產生的機械碼來說: ([Pentium 4](https://en.wikipedia.org/wiki/Pentium_4)) > 資料來源: [The Software Optimization Cookbook: High Performance Recipes for IA-32 Platforms, 2nd Edition](http://www.amazon.com/The-Software-Optimization-Cookbook-Performance/dp/0976483211) > 針對 2004 年的個人電腦 |指令|數量|延遲|cycle count| |:---:|:---:|:---:|:----:| |減法|5 |1 |5 | |比較|5 |1 |5 | |分支|14 |1 |14 | |其他|0 |1 |0 | | | | |24 | 改成 mod 操作的 C 語言程式如下： (v2) ```cpp int findGCD(int a, int b) { while (1) { a %= b; if (a == 0) return b; if (a == 1) return 1; b %= a; if (b == 0) return a; if (b == 1) return 1; } } ``` 以 Intel C compiler 產生的機械碼來說: |指令 |數量|延遲|cycle count| |:-----------:|:--:|:--:|:--------:| |mod (整數除法) |2 |68 |136 | |比較 |3 |1 |3 | |分支 |3 |1 |3 | |其他 |6 |1 |6 | | | | |148 | 若可排除整數除法和減法的負擔，得到以下實做： (v3) ```cpp int findGCD(int a, int b) { while (1) { if (a > (b * 4)) { a %= b; if (a == 0) return b; if (a == 1) return 1; } else if (a >= b) { a -= b; if (a == 0) return b; if (a == 1) return 1; } if (b > (a * 4)) { b %= a; if (b == 0) return a; if (b == 1) return 1; } else if (b >= a) { b -= a; if (b == 0) return a; if (b == 1) return 1; } } } ``` 可得到更快的實做: ``` v1: 14.56s v2: 18.55s v3: 12.14s ``` 上述都是針對 2004 年的個人電腦所得的數據。 ### 針對 [Arm Cortex-A8](https://en.wikipedia.org/wiki/ARM_Cortex-A8) 進行實驗 * 平台：BeagleBone Black (ARMv7-A, 1 GHz) * 實驗程式碼：[GCD_Efficiency_Compare](https://github.com/JaredCJR/GCD_Efficiency_Compare) 採用組合語言撰寫 (GCD 的部份)，由於 Cortex-A8 不支援 idiv (要在 ARM Cortex-A15 或 A7 以上才有)，因此改用以 ARM 組合語言實做 * 發現 gcc 的除法會針對 2 的倍數做判斷，移位（除法） * 實做的版本尚未考慮 exception - [A Fast Hi Precision Fixed Point Divide](http://me.henri.net/fp-div.html) 列出 exception 重點 ### 實驗方法 * 圖上一個點，代表 GCD (big, small) 中的 big，會尋找比它小的所有正整數（>1）的GCD所需的時間「總和」 * 若存在 big=9995，那麼會找 9993 組 GCD，small 的範圍從 2 到 9994 * 如: (9995,2), (9995,3), (9995,4),......, (9995,9994) * **所以 big 越大，找尋的組數越多，整體趨勢微微上升（從左到右）是正常的** * 重點是比較各演算法的趨勢 ### ARM vs.  Thumb instruction * ARM instruction速度較快 * 直線為 curve fitting 的結果(時間越低，效能越高) ![](https://hackpad-attachments.s3.amazonaws.com/embedded2015.hackpad.com_PxVP6UWvGxh_p.339782_1439772172448_undefined) **Assembly** ```cpp MOD: @ Entry r0: low (remainder low) must be postive<input:Dividend> @ r1: high (remainder high) any number @ r2: Divisor (divisor) must be non-zero and @ positive<input:Divisor> low .req r0; high .req r1; divisor .req r2; mov high,#0 .rept 33 @ repeat 33 times subs high, high, divisor @ remainder = remainder - divisor it lo addlo high, high, divisor @ low << 1; ((high-divisor) > 0) ? (new bit = 1) : (new bit = 0) adcs low, low, low @ divisor >> 1 == high(64bit for high & low) << 1 adc high, high, high .endr @ shift remainder to correct position @ because last divide don't know this is the last divide, @ it still prepare for the next divide<shift>) lsr high, #1 @ Exit r0: Quotient (remainder low) @ r1: remainder (remainder high) mov pc, lr ``` **Euclidean(thumb) V.S. mod(thumb) V.S. mod_minus(thumb)** 不符合預期的原因推估 * v3 版本的 branch 過多（if,else），thumb 的 branch 效能不佳，導致拉低平均 * if 數量比較（一個 loop） - v1 : 3 個 - v3 : 7 個 ![](https://hackpad-attachments.s3.amazonaws.com/embedded2015.hackpad.com_PxVP6UWvGxh_p.339782_1439770375774_undefined) **Euclidean(arm) V.S. mod(arm) V.S. mod_minus(arm)** 比較 * 與 thumb 版本 * arm 的 branch 效能較佳，所以 v3 版本速度改善最多，這裡符合預期 * v1 輸給 v2 * arm 的 v1 與 thumb 的 v1 相比，只有微微速度提昇（看第一張圖），但是 mod 卻提昇較多（幾乎是半張圖片的 range），導致這個結果 * curve fitting ![](https://hackpad-attachments.s3.amazonaws.com/embedded2015.hackpad.com_PxVP6UWvGxh_p.339782_1439741803792_undefined) ## C 語言 * arm-angstrom-linux-gnueabi-gcc 內的除法經過最佳化後，效能比減法更佳 ![](https://hackpad-attachments.s3.amazonaws.com/embedded2015.hackpad.com_PxVP6UWvGxh_p.339782_1439777592642_embedded2015.hackpad.com_PxVP6UWvGxh_p%20(1).png) * arm-linux-gnueabihf-gcc(linaro) (without vfpv3) * 上下圖比較，可看到整體在 arm-angstrom-linux-gnueabi-gcc 效能比較好 * 即使在這裡，單純除法比單純減法（的演算法）還快... * 可是綜合起來看 v3 更快(符合預測) ![](https://hackpad-attachments.s3.amazonaws.com/embedded2015.hackpad.com_PxVP6UWvGxh_p.339782_1440082584221_undefined) * arm-linux-gnueabihf-gcc(linaro) (with vfpv3) * mod 加快，減法變慢 * 開了導致反效果 ![](https://hackpad-attachments.s3.amazonaws.com/embedded2015.hackpad.com_PxVP6UWvGxh_p.339782_1440118082155_undefined) ## [Binary GCD](https://en.wikipedia.org/wiki/Binary_GCD_algorithm) > [Binary GCD Algorithm](https://iq.opengenus.org/binary-gcd-algorithm/) binary gcd 的精神就是當兩數為偶數時，必有一 $2$ 因子。 $$\text{gcd}(x, y) = 2·\text{gcd}(\frac{x}{2}, \frac{y}{2})$$ 且一數為奇數另一數為偶數，則無 $2$ 因子。 $$ \text{gcd}(x, y) = \text{gcd}(x, \frac{y}{2})$$ 於是我們可以改良為： $$ \text{even_factor} = \text{min}(\text{ctz}(x), \text{ctz}(y)) \\ \text{gcd}(x, y) = \text{even_factor}·\text{gcd}((x\gg \text{even_factor}), (y\gg \text{even_factor})) $$ 其中符號 $\gg$ 是 right shift。 剩下的過程就直接採用 Euclidean algorithm 的作法即可。 注意，僅在微處理器具備有效的 `ctz` 指令時，[Binary GCD](https://en.wikipedia.org/wiki/Binary_GCD_algorithm) 才會比 Euclidean 來得快。 延伸閱讀: * [The Speed of GCD](https://hbfs.wordpress.com/2013/12/10/the-speed-of-gcd/) * [Fastest way to compute the greatest common divisor](https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/) ## Optimized Binary GCD for Modular Inversion > it requires GCC or Clang, a 64-bit x86 architecture, and a target CPU that supports the [BMI2 opcodes](https://en.wikipedia.org/wiki/Bit_manipulation_instruction_set) (i.e. Haswell or later, in the line of Intel CPU). [bingcd](https://github.com/pornin/bingcd) [^factor]: 整數 $n$ 的**因數**為可**整除** $n$ 的數 [^euclidean]: [Wikipedia/ 輾轉相除法的展示動畫](https://zh.wikipedia.org/zh-tw/%E8%BC%BE%E8%BD%89%E7%9B%B8%E9%99%A4%E6%B3%95#/media/File:Euclidean_algorithm_252_105_animation_flipped.gif)
Published on HackMD
1696