最大公約数を求める2つのアルゴリズムを書いてみた

Cで多倍長整数演算のプログラムを書くのに熱中していた。既存のライブラリ*1を参考にしながら、自分でコードを書いてみている。勉強になるなあ。

アルゴリズムの参考にしたのは、クヌース先生の『The Art of Computer Programming Volume 2 Seminumerical Algorithms』。バイブルの名は伊達ではないようで。どれだけの知見を持っていたらこんなものが書けるのか……読み物としても結構面白い。

ただ、まだ乗算とか面倒な部分は実装していないので、本格的に役に立つのはまだ先。Karatsuba法とかFFTを使った高速乗算法、うまく書けるかなあ。

互除法ともう一つ

ここで話が変わる。TAOCPの他の部分を読んでいたら、最大公約数を求めるアルゴリズムが載っていた。有名なユークリッドの互除法の他に、もう一つ。今ここに書き残しておかないと忘れる自信があるので、忘れる前に書いておく。

まず、ユークリッドの互除法。

int gcd_euclid(int u, int v) {
  int r;
  while (0 != v) {
    r = u % v; u = v; v = r; /* swap */
  }
  return u;
}

きれいだ。この方法は十分に簡潔で高速なように思われる。

ところが、まだあるらしい。

Euclidによる始祖的アルゴリズムは何世紀にもわたって使われてきているから、実はそれが必ずしも最大公約数を求める最良の方法ではないという事実にはかなり驚く。
主に二進の算術演算に適した全く別のアルゴリズムをJosef Steinが1961年に考案している。

これもCで書いてみた。

int gcd_stein(int u, int v) {
  int k, t;

  /* u, vどちらかが奇数になるまで2で割り続ける */
  k = 0;
  while((0 == (u & 1)) && (0 == (v & 1))) {
    k++;
    u /= 2;
    v /= 2;
  }

  if ((u & 1) == 0) {
    t = u / 2;
  } else {
    t = -v;
  }

  do {
    /* tが偶数である限り、2で割り続ける */
    while ((t & 1) == 0) t /= 2;

    if (t > 0) {
      u = t;
    } else {
      v = -t;
    }

    t = u - v;
  } while (0 != t);
  return u * (1 << k);
}