拡張ユークリッドの互除法を理解する

拡張ユークリッドの互除法とは？

${\rm gcd}(a,b)$ を求めるついでに $xa + yb = {\rm gcd}(a,b)$ の解 $(x, y)$ もついでに求めるアルゴリズムです

ユークリッドの互除法復習

とりあえず147と60のgcdを求めてみましょう

\begin{alignat}{2}
& & & & 147 \\
& & & & 60 \\
147 & \quad \% \quad & 60 \quad = \quad & & 27 \\
60 & \quad \% \quad & 27 \quad = \quad & & 6 \\
27 & \quad \% \quad & 6 \quad = \quad & & 3 \\
6 & \quad \% \quad & 3 \quad = \quad & & 0 \\
\end{alignat}

よって ${\rm gcd}(147,60) = 3$ です
見ての通り $z_0 = 147, z_1 = 60, z_i = z_{i-2} \% z_{i-1}$ という隣接2項間漸化式ですね
フィボナッチ数列みたいな感じです

これを素直に実装するとすれば以下のような感じでしょうか

int gcd(int a, int b){
    int z[1000] = { a, b };
    for(int i = 2; ; ++i){
        int q = z[i-2] / z[i-1];
        z[i] = z[i-2] - z[i-1] * q;
        if(z[i] == 0) return z[i-1];
    }
}

後の説明が楽になるようにmod演算ではなく除算と引き算を使って実装しました

拡張ユークリッドの互除法

ここで $z_i = z_{i-2} - z_{i-1} * q_i, \quad q_i = {\rm floor}(z_{i-2} / z_{i-1})$ という処理を
$xa + yb = z$ という形の数式に対して実行してみましょう

\begin{alignat}{5}
1 & \quad \times \quad & 147 & \quad + \quad & 0 & \quad \times \quad & 60 & \quad = \quad & 147 & \qquad (1) \\
0 & \quad \times \quad & 147 & \quad + \quad & 1 & \quad \times \quad & 60 & \quad = \quad & 60 & \qquad (2) \\
1 & \quad \times \quad & 147 & \quad + \quad & (-2) & \quad \times \quad & 60 & \quad = \quad & 27 & \qquad (3) = (1) - 2 \times (2) \\
(-2) & \quad \times \quad & 147 & \quad + \quad & 5 & \quad \times \quad & 60 & \quad = \quad & 6 & \qquad (4) = (2) - 2 \times (3) \\
9 & \quad \times \quad & 147 & \quad + \quad & (-22) & \quad \times \quad & 60 & \quad = \quad & 3 & \qquad (5) = (3) - 4 \times (4) \\
& & & \quad \text{(略)} &  &  &  & \quad = \quad & 0 & \qquad (6) = (7) - 2 \times (8) \\
\end{alignat}

ということで
${\rm gcd}(147,60) = 3$ を求めるついでに $xa + yb = {\rm gcd}(147,60) = 3$ の解 $(x, y) = (9, -22)$ が求まりました。

先ほどと同じように実装するとこんな感じでしょうか？

struct eq {
    int x,y,z;
    eq operator*(int q)const{ return { x*q, y*q, z*q }; }
    eq operator-(const eq &b)const{ return { x-b.x, y-b.y, z-b.z }; }
};
eq extgcd(int a, int b){
    eq e[1000] = { {1,0,a}, {0,1,b} };
    for(int i = 2; ; ++i){
        int q = e[i-2].z / e[i-1].z;
        e[i] = e[i-2] - e[i-1] * q;
        if(e[i].z == 0) return e[i-1];
    }
}

フィボナッチ数列を求める時と同じように過去2回分まで結果を持っておけば十分なので整理すると

std::tuple<int,int,int> extgcd(int a, int b){
    int x=1,y=0;
    int u=0,v=1;
    while(b){
        int q = a / b;
        x -= u * q;
        y -= v * q;
        a %= b;
        std::swap(x, u);
        std::swap(y, v);
        std::swap(a, b);
    }
    return {x,y,a};
}