拡張ユークリッド互除法の非再帰的実装を理解する(C++)

この記事は僕が拡張ユークリッド互除法の非再帰的実装を理解した道筋をメモしたもの。

拡張ユークリッド互除法について

$\displaystyle{ ax + by = \mathrm{gcd}(a,b) }$

の整数解 (x, y) をひとつ求めることができるアルゴリズム。高校の数学Aで初めて習ったけれど、こんな名前がついてることは最近知った。

拡張ユークリッド互除法を a = 493、b = 408 のときの具体例を使って簡単に説明する。まず、493、408に対してユークリッド互除法を行う。

$\displaystyle{ \begin{aligned} 493 &= 1 \cdot 408 + 85\\ 408 &= 4 \cdot 85 + 68\\ 85 &= 1 \cdot 68 +17\\ 68 &= 4 \cdot 17 + 0 \end{aligned} }$

これで、gcd(493,408) = 17 が分かった。上の一連の結果を利用すると

$\displaystyle{ \begin{aligned} 17 &= 85 - 1 \cdot 68\\ &= (493 - 1 \cdot 408) - 1 \cdot (408 - 4 \cdot 85)\\ &= (493 - 1 \cdot 408) - 1 \cdot \{ 408 - 4 \cdot (493 - 1 \cdot 408) \}\\ &= 493 \cdot 5 + 408 \cdot (-6) \end{aligned} }$

が分かり、整数解 (x, y) = (5, -6) が得られる。

これをC++で非再帰的に実装することを目指す。でも、実際に実装しようとするとう～ん？となる。しくみは単純なのに...... 頭のいい人なら「こうすればいいじゃん！」ってすぐになるのかな？ぼくはならない。

拡張ユークリッド互除法を式で表現する

とりあえず、拡張ユークリッド互除法がやっていることを式で表現してから考えようと思った。そして、Wikipediaを見るとその答えが書いてあった。いろいろサイトを見て回ったが、個人的にはWikipediaが一番わかりやすかった。ほぼWikipediaのままになるが、拡張ユークリッド互除法の過程を式で表現していく。

まず始めにユークリッド互除法を行う。 $r_0 = a,\ r_1 = b$ として次のようにあまりが0になるまで割り算を繰り返す。

$\displaystyle{ \begin{aligned} r_0 &= k_0 r_1 + r_2\\ r_1 &= k_1 r_2 + r_3\\ &\ \vdots\\ r_{h-2} &= k_{h-2} r_{h-1} + r_{h}\\ r_{h-1} &= k_{h-1} r_{h} + 0 \end{aligned} }$

ただし、 $\displaystyle k_i = \left\lfloor \frac{r_{i}}{r_{i+1}} \right\rfloor$ 。上の式は次のように書き直すことができる。

$\displaystyle{ \begin{aligned} \begin{bmatrix} r_0 \\ r_1 \end{bmatrix} &= \begin{bmatrix} k_0 & 1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} r_1 \\ r_2 \end{bmatrix}\\ \begin{bmatrix} r_1 \\ r_2 \end{bmatrix} &= \begin{bmatrix} k_1 & 1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} r_2 \\ r_3 \end{bmatrix}\\ &\ \vdots\\ \begin{bmatrix} r_{h-1} \\ r_{h} \end{bmatrix} &= \begin{bmatrix} k_{h-1} & 1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} r_{h} \\ 0 \end{bmatrix} \end{aligned} }$

そして、これらの式を連結することで次の式が得られる。

$\displaystyle{ \begin{aligned} \begin{bmatrix} r_0 \\ r_1 \end{bmatrix} = \begin{bmatrix} k_0 & 1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} k_1 & 1\\ 1 & 0 \end{bmatrix} \cdots \begin{bmatrix} k_{h-1} & 1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} r_{h} \\ 0 \end{bmatrix} \end{aligned} }$

ここで、

$\displaystyle{ \begin{bmatrix} k_{i} & 1\\ 1 & 0 \end{bmatrix}^{-1} = \begin{bmatrix} 0 & 1\\ 1 & -k_{i} \end{bmatrix} }$

と $r_0 = a,\ r_1 = b,\ r_{h} = \mathrm{gcd}(a,b)$ を考慮すると次の式が得られる。

$\displaystyle{ \begin{aligned} \begin{bmatrix} \mathrm{gcd}(a,b) \\ 0 \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & -k_{h-1} \end{bmatrix} \begin{bmatrix} 0 & 1 \\ 1 & -k_{h-2} \end{bmatrix} \cdots \begin{bmatrix} 0 & 1 \\ 1 & -k_{0} \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} \end{aligned} }$

上の式より、

$\displaystyle{ \begin{aligned} \begin{bmatrix} x & y \\ u & v \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & -k_{h-1} \end{bmatrix} \begin{bmatrix} 0 & 1 \\ 1 & -k_{h-2} \end{bmatrix} \cdots \begin{bmatrix} 0 & 1 \\ 1 & -k_{0} \end{bmatrix} \end{aligned} }$

とおけば、 $(x, y)$ は $ax + by = \mathrm{gcd}(a,b)$ を満たし、解が求まった。

実装する

拡張ユークリッド互除法の過程を式で表現したことによって、実装へのハードルがかなり下がった。上の式を眺めれば、とりあえず、次のようにすれば拡張ユークリッド互除法を行えることがわかる。

Step0.(初期化)
$i \leftarrow 0,\ r_0 \leftarrow a,\ r_1 \leftarrow b$

$\displaystyle{ \begin{bmatrix} x & y \\ u & v \end{bmatrix} \leftarrow \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} }$

Step1.
$k_i \leftarrow \displaystyle \left \lfloor \frac{r_{i}}{r_{i+1}} \right \rfloor$
$r_{i+2} \leftarrow r_{i} - k_{i} r_{i+1}$

Step2.

$\displaystyle{ \begin{bmatrix} x & y \\ u & v \end{bmatrix} \leftarrow \begin{bmatrix} 0 & 1 \\ 1 & -k_{i} \end{bmatrix} \begin{bmatrix} x & y \\ u & v \end{bmatrix} }$

Step3.
もし $r_{i + 2} = 0$ ならば現在の $x, y$ を解として終了。
そうでなければ $i \leftarrow i + 1$ して Step.1へ移動。

さて、上の手順をもうすこし整理する。上では $\{k_{i}\}$ や $\{r_{i}\}$ といった配列を使っているが、 $\{k_{i}\}$ は $k$ という1つの変数で大丈夫だし、 $\{r_{i}\}$ もうまくすれば2つの変数を使うだけでなんとかなりそうだ。実際、 $r_{i}$ は $k_{i}$ の計算が終わるともう使わないから、その領域に $r_{i+2}$ を上書きすれば良い。Step2. の行列計算は書き下すと

$\displaystyle{ \begin{alignat}{2} x &\leftarrow u, &\quad y &\leftarrow v, \\ u &\leftarrow x - k_i u, &\quad v &\leftarrow y - k_i v \end{alignat} }$

となっている。以上を踏まえると、アルゴリズムは次のように整理できる。

Step0.(初期化)
$s \leftarrow a,\ t \leftarrow b$
$x \leftarrow 1,\ y \leftarrow 0,\ u \leftarrow 0,\ v \leftarrow 1$

Step1.
$k \leftarrow \displaystyle \left \lfloor \frac{s}{t} \right \rfloor$
$s \leftarrow s - kt$
$\mathrm{swap}(s,t)$

Step2.
$x \leftarrow x - ku$
$\mathrm{swap}(x,u)$
$y \leftarrow y - kv$
$\mathrm{swap}(y,v)$

Step3.
もし $t = 0$ ならば現在の $x, y$ を解として終了。
そうでなければStep.1へ移動。

実装例(C++)

template <typename T>
constexpr std::pair<T, T> extended_euclidean(const T& a, const T& b) {
    T x = 1, y = 0, u = 0, v = 1, s = a, t = b;
    T k = 0, tmp = 0;
    while (t) {
        k = s / t, s -= k * t;
        tmp = s, s = t, t = tmp;
        x -= k * u;
        tmp = x, x = u, u = tmp;
        y -= k * v;
        tmp = y, y = v, v = tmp;
    }
    return {x, y};
}

モジュラ逆数を求める(応用)

拡張ユークリッド互除法はモジュラ逆数を求めるのにも使われる。フェルマーの小定理を使う方法と違って法は素数でなくてもいいし、さらに速いらしい。詳細はWikipediaなどを参照。モジュラ逆数を求めるときは $y$ を求めなくていいから、 $y, v$ についての計算を行う必要がない。そのため、実装は次のようになる。

template <int64_t Modulus>
constexpr int64_t modinv(const int64_t& a) {
    int64_t x = 1, u = 0, s = a, t = Modulus;
    int64_t k = 0, tmp = 0;
    while (t) {
        k = s / t, s -= k * t;
        tmp = s, s = t, t = tmp;
        x -= k * u;
        tmp = x, x = u, u = tmp;
    }
    x %= Modulus;
    if (x < 0) x += Modulus;
    return x;
}

スコルの知恵袋

主にプログラミング関係の気になったこと等をまとめたブログ