[程式] 不用除法的最大公因數 (GCD)

有不少數論 (Number Theory) 的計算中都會用到找兩個正整數的最大公因數 (greatest common divisor, GCD)，而最基本常用找最大公因數的方法中就是利用輾轉相除法 (Euclidean algorithm)，如下

$gcd(x, y) = gcd(y, x - q \times y)$

其中 q 為任意整數，下面的這段程式實作出輾轉相除法的概念：

int gcd(int x, int y)
{
 while ((x %= y) && (y %= x)) ;
 return (x + y);
}

然而，在現在處理器中求餘數的運算時間相較於基本的加減法以及位元運算 (bitwise operation) 是慢很多的，因此 Knuth 的 TAOCP (The Art of Computer Programming) 中的第 4.5.2 節就有提到其實曾經有人提出了另一種方法，僅僅使用位元運算以及加減法來找最大公因數。

這個方法的基本概念如下：

• 若 x 跟 y 都是偶數，則：$gcd(x, y) = 2 \times gcd(x/2, y/2)$
• 若 x 是偶數而 y 是奇數，則：$gcd(x, y) = gcd(x / 2, y)$
• 若 x 跟 y 都是奇數，則：$x - y$ 是偶數，且 $|x - y| < max(x, y)$
• 同輾轉相除法，不過簡化成：$gcd(x, y) = gcd(x - y, y)$

透過以上四點，其實作如下：

int gcd(int x, int y)
{
 int k = 0;
 for (; ((x & 1) == 0) && ((y & 1) == 0); ++k) { x >>= 1; y >>= 1; }
 
 int t = (x & 1) ? (-y) : x;
 while (t) {
  while ((t & 1) == 0) { t >>= 1; }
  if (t > 0) { x = t; }
  else { y = -t; }
  t = x - y;
 }
 return (x << k);
}

• 第 1 ~ 2 行是為了先讓至少其中一個數字變為奇數，同時我們也能知道 $2^k$ 是 x 與 y 的公因數，但 $2^{k+1}$ 就不是 x 與 y 的公因數
• 接下來直到 return 前則是類似輾轉相除法的過程，持續運作到 x 跟 y 相等，最後的最大公因數即為 $2^k \times x$。
• 當 x 與 y 不同時，利用輾轉相除法將 x、y 中較大者重設為 $x - y$
• 當 x 或 y 其中一個是偶數時，因為 2 必定不是 x 或 y 的公因數，因此可以一直除以 2 直到變為奇數。

以下的這個連結分別實作出兩種方法並且比較效能：
https://ideone.com/nhu1HV

不過就結果來說，可以發現反而還比較慢...這大概是因為現在的處理器對於除法的運算其實並不慢，相較於第二種版本，雖然只使用了 bitwise operation 跟減法，然而它所需的迴圈數量、指令數目甚至於條件判斷都多於輾轉相除法，因此除非除法運算的速度慢到難以忍受，不然直接實作出輾轉相除法其實還比較簡單。