Try   HackMD

最大公因數特性和實作考量

貢獻者: jserv, TibaChang, Tcc0403

簡介

本文探討最大公因數 (Greatest Common Divisor, GCD) 特性，考慮微處理器架構實作帶來的效能衝擊，並探討如何運用 Linux 核心提供的 perf 工具進行效能分析並改進。

Greatest Common Divisor 中譯「最大公因數」，簡稱 GCD。
給定整數 \(a, b\)，二者各自有自己的因數[1]，取相同的因數中最大的數，即最大公因數 \(\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

歐幾里得法，也就是「輾轉相除法」

給定整數 \(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} \]

[2]

展示動畫

綜合上面過程，可得程式碼實作：

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

給定整數 \(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} \]

綜合上述過程：

/* 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;
}

延伸閱讀:

GCD 實作和微處理器指令

回顧歐幾里得法:

  1. 大數 := 大數 - 小數
  2. 若「最大」和「最小」兩數相同，那麼它就是最大公因數，否則，回到第 1 步

示範 GCD(48, 40) 的計算過程

  • 48 - 40, 40 \(\to\) 8, 40
  • 8 != 40, 於是重複上述動作: 40 - 8, 8 \(\to\) 32, 8
  • 8 != 32, 於是重複上述動作: 32 - 8, 8 \(\to\) 24, 8
  • 8 != 24, 於是重複上述動作: 24 - 8, 8 \(\to\) 16, 8
  • 8 != 16, 於是重複上述動作: 16 - 8, 8 \(\to\) 8, 8
  • 8 == 8 成立，於是 8 為 GCD

用 C 語言實現歐幾里得法: (v1)

int findGCD(int a, int b)
{
    while (1) {
        if (a > b) a -= b;
        else if (a < b) b -= a;
        return a;
    }
}

針對 Pentium 4 硬體的實驗

以 Intel C compiler 產生的機械碼來說: (Pentium 4)

資料來源: The Software Optimization Cookbook: High Performance Recipes for IA-32 Platforms, 2nd Edition
針對 2004 年的個人電腦

指令 數量 延遲 cycle count
減法 5 1 5
比較 5 1 5
分支 14 1 14
其他 0 1 0
24

改成 mod 操作的 C 語言程式如下： (v2)

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)

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 進行實驗

採用組合語言撰寫 (GCD 的部份)，由於 Cortex-A8 不支援 idiv (要在 ARM Cortex-A15 或 A7 以上才有)，因此改用以 Arm 組合語言實作

實驗方法

  • 圖上一個點，代表 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 指令執行較快
  • 直線為 curve fitting 的結果 (時間越低，效能越高)

image

Assembly

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 個

image

Euclidean(arm) V.S. mod(arm) V.S. mod_minus(arm)

比較

  • 與 thumb 版本
    • arm 的 branch 效能較佳，所以 v3 版本速度改善最多，這裡符合預期
  • v1 輸給 v2
    • arm 的 v1 與 thumb 的 v1 相比，只有些許提昇 (對照第一張圖），但是 mod 卻提昇較多（幾乎是半張圖片的範圍），導致這個結果
  • curve fitting
    image

C 語言

  • arm-angstrom-linux-gnueabi-gcc 內的除法經過最佳化後，效能比減法更佳

image

  • arm-linux-gnueabihf-gcc(linaro) (without vfpv3)
    • 上下圖比較，可看到整體在 arm-angstrom-linux-gnueabi-gcc 效能比較好
      • 即使在這裡，單純除法比單純減法（的演算法）還快…
      • 可是綜合起來看 v3 更快(符合預測)

image

  • arm-linux-gnueabihf-gcc(linaro) (with vfpv3)
    • mod 加快，減法變慢
    • 開了導致反效果

image

Binary GCD

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 才會展現相較於 Euclidean 的優勢，這也是為何 Linux 核心提供二種 gcd 實作。

unsigned long gcd(unsigned long a, unsigned long b)
{
    unsigned long r = a | b;
    if (!a || !b) return r;
    b >>= __builtin_ctzl(b);
    if (b == 1) return r & -r;
    for (;; a -= b) {
        a >>= __builtin_ctzl(a);
        if (a == 1) return r & -r;
        if (a == b) return a << __builtin_ctzl(r);
        if (a < b) { unsigned long t = a; a = b; b = t; }
    }
}

延伸閱讀:

改進實作效能

bingcd

it requires GCC or Clang, a 64-bit x86 architecture, and a target CPU that supports the BMI2 opcodes (i.e. Haswell or later, in the line of Intel CPU).

考慮以下實作：

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int shift;
    for (shift = 0; !((u | v) & 1); shift++) {
        u /= 2, v /= 2;
    }
    while (!(u & 1))
        u /= 2;
    do {
        while (!(v & 1))
            v /= 2;
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v);
    return u<<shift;
}

GCD 演算法

  1. If both x and y are 0, gcd is zero \(gcd(0, 0) = 0\).
  2. \(gcd(x, 0) = x\) and \(gcd(0, y) = y\) because everything divides 0.
  3. If x and y are both even, \(gcd(x, y)=2 \times gcd(\dfrac{x}{2},\dfrac{y}{2})\) because \(2\) is a common divisor. Mulitplication with \(2\) can be done with betwise shift operator.
  4. If x is even and y is odd, \(gcd(x, y) = gcd(\dfrac{x}{2}, y)\).
  5. On similar lines, if x is odd and y is even, then \(gcd(x, y) = gcd(x, \dfrac{y}{2})\). It is because \(2\) is not a common divisor.

以下分段說明：

if (!u || !v) return u | v;

\(\to\) 對應至 GCD 演算法第 2 點

  • \(gcd(x,0)=x\) and \(gcd(0, y)=y\) because everything divides 0.

若 u 和 v 中其中一數為 0 ，則直接回傳另一數

x | 0 會得到 x

int shift;
for (shift = 0; !((u | v) & 1); shift++) {
    u /= 2, v /= 2;
}

\(\to\) 對應至 GCD 演算法第 3 點

  • If x and y are both even, \(gcd(x, y) = 2 \times gcd(\dfrac{x}{2},\dfrac{y}{2})\) because \(2\) is a common divisor. Mulitplication with \(2\) can be done with betwise shift operator.

先取出u 和 v 之間屬於 \(2^n\) 的公因數，shift 即為最大的 \(n\)

while (!(u & 1))
    u /= 2;

\(\to\) 對應至 GCD 演算法第4 點

  • If x is even and y is odd, \(gcd(x, y) = gcd(\dfrac{x}{2}, y)\).

當 u 為偶數時，可以除以 \(2\)

do {
    while (!(v & 1))
        v /= 2;
    if (u < v) {
        v -= u;
    } else {
        uint64_t t = u - v;
        u = v;
        v = t;
    }
} while (v);

\(\to\) 其中 while(!(v &1)) v /= 2; 對應至 GCD 演算法第 5 點

  • On similar lines, if x is odd and y is even, then \(gcd(x, y) = gcd(x, \dfrac{y}{2})\). It is because \(2\) is not a common divisor.

可以透過迴圈不變性 (loop invariant) 來幫助解讀程式

  • u 為被除數
  • v 為除數

當除數 v為 0 時結束迴圈，被除數 u 為進入迴圈前 u 和 v 的最大公因數

事實上 u 在迴圈前不一定是被除數，因為先前透過 GCD 演算法第 4 點去除不必要的計算

return u << shift;

最後回傳最大公因數時，乘上原先取出屬於 \(2^n\) 的公因數，可以透過左移達成

透過 __builtin_ctz 改寫程式

已知除法運算會耗費較多時間，我們可以利用 __builtin_ctz (編譯器會產生對應到現代微處理器的專門指令) 搭配上述針對二進位特性調整的演算法，減少除法的使用

將一百萬組由 rand 函式生成的亂數傳入 gcd64 函式執行，並藉由 Linux 核心提供的 perf 工具來分析：

$ gcc -g3 -O0 gcd.c -o gcd
$ perf record ./gcd
$ perf annotate

image

由上圖可以發現，花費 CPU 週期數最多的段落為迴圈當中的 v /= 2 運算，佔據約 28%

以下為改寫過後的程式

uint64_t gcd64_ctz(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int shift = __builtin_ctzll(u | v);
    u = u >> __builtin_ctzll(u);
    
    do {
        v = v >> __builtin_ctzll(v);
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v);
    return u << shift;
}

透過 perf 分析前後差異

$ gcc -g -O0 gcdcmp.c -o gcdcmp
$ perf record ./gcdcmp
$ perf report

image

將同樣一百萬筆由 rand 函式所產生的亂數傳入 gcd64 和 gcd64_ctz 函式比較，發現改寫過後的版本花費時間幾乎是原來的一半。

透過 perf 分析改寫程式

$ gcc -g -O0 gcd_ctz.c -o gcd_ctz
$ perf record ./gcd_ctz
$ perf annotate

image

佔據花費時間最多的為分支指令 (如 jne)，其次是迴圈內部的減法運算

  1. 整數 \(n\) 的因數為可整除 \(n\) 的數 ↩︎

  2. Wikipedia/ 輾轉相除法的展示動畫 ↩︎

Last changed by 
 
sysprog
System Software Programming
1
6948

Read more

你所不知道的 C 語言：數值系統

回歸第一手資料，透過反思 C 語言程式設計的細節，重新學習電腦原理

Jul 21, 2025
建構 User-Mode Linux 的實驗環境

自 Linux 核心原始程式碼編譯出 User-mode Linux 並運用工具來建構實驗所需的檔案系統

Jul 20, 2025
Linux 核心設計: Scheduler(8): Energy Aware Scheduling

系統的 throughput 是評斷排程器優劣的重要效能

Jul 19, 2025
你所不知道的 C 語言: linked list 和非連續記憶體

無論是作業系統核心、C 語言函式庫內部、程式開發框架，到應用程式，都不難見到 linked list 的身影，包含多種針對效能和安全議題所做的 linked list 變形，又還要考慮到應用程式的泛用性 (generic programming)，是很好的進階題材。

Jul 17, 2025
Read more from sysprog
Published on HackMD