「フゥーハハハ！　これが運命石の扉の選択か……
*　これより、オペレーション・GCDを開始する！！」*

前回、C++11で黒魔術的NTZとNLZを実装する話をやりました。今回はその続きです。
まあ、続きと言っても、NTZを使えばこういうことができるよ！　ってだけなのですが。

SteinのGCDアルゴリズム

タイトルと冒頭で謎のシュタゲネタをぶち込みましたが、今回実装するのは「バイナリGCDアルゴリズム」とか「Steinのアルゴリズム」と呼ばれるアルゴリズムです。
Steinというのは別に運命石ではなくて、人名です。Josef Steinさんという方が発見した最小公倍数アルゴリズムらしいです。
前回と同じくネタ元のブログ記事を紹介します。

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

前回のNTZ/NLZと比較するとそれほど黒魔術的ではありません。
考え方はこうです。

ユークリッドの互除法は美しいが最速ではない。
なぜ速くないか→剰余が遅いから
なら剰余をを使わなければいい
剰余の代わりに全部シフト演算にできたら速くなるんじゃね？

上記のブログ記事では非常に素直な実装を行っていますが、私のような狂気のマッドC++erはこう考えます。

ループ回してシフト演算してる所、NTZに置き換えられるんじゃね？
C++11らしくconstexprにしたいよね。再帰しかないな！

やってみました。

gcd_stein.hpp

# include <ntz_nlz.hpp>

// std::minがconstexprじゃないのが悪い
template<typename T> inline constexpr T min(T x, T y) noexcept { return x<y?x:y; }

template<typename T>
struct gcd_stein_
{
    using result_type = decltype(std::declval<T>()%std::declval<T>());
    static inline constexpr result_type func1(T x, T y, int n, int m) noexcept { 
        return func2(x>>n, y>>m) << min(n, m);
    }
    static inline constexpr result_type func2(T x, T y) noexcept {
        return x==y? x: (x>y ? func2(func3(x-y),y) : func2(func3(y-x),x));
    }
    static inline constexpr result_type  func3(T x) noexcept { return x >> ntz(x); }
};

template<typename T>
inline constexpr typename gcd_stein_<T>::result_type gcd_stein(T x, T y) noexcept {
     return gcd_stein_<T>::func1(x, y, ntz(x), ntz(y));
}

意外と短い。そして元記事のサンプルコードが影も形もない。
なお、ここで実装したGCDは、x, yのどちらかが0以下であった場合を考慮していません。あくまで正の2つの整数の最大公約数を求める関数であることをご了承ください。

それでは、短いですが本日はここまで。