拡張ユークリッドの互除法の仕組みを非再帰版で理解する

整数の問題で方程式 $a x + b y = d$ の解を求めたりモジュラ逆数を計算したりする際、「拡張ユークリッドの互除法」 (extended Euclidean algorithm) によくお世話になります。ただ、このアルゴリズムがどういう仕組みなのか、いくつか解説を見ても難しくて飲み込めずにいました。

そんな中、 Wikipedia 英語版の記事（特に擬似コード）を見たところ、「拡張」のしかたがとても単純で一気に理解が進みました。本記事では実際に手計算して仕組みを理解し、ついでに Python で実装してみます。

なお、本記事ではアルゴリズムの仕組みを示すだけとし、正しさの証明などについては触れません。

TL;DR

方程式 $a x + b y = \gcd (a, b)$ の整数解 $x, y$ を求めたいときは、いきなり完全な解を目指すのではなく、 $a x_k + b y_k = d_k$ を満たす $x_k, y_k, d_k$ を見つけて改善していきます。

以下の2つの自明な式を出発点に、右辺が通常の互除法になるよう式同士で引き算を繰り返すと答えが現れます。

\begin{align}
a \cdot 1 + b \cdot 0 &= a \\
a \cdot 0 + b \cdot 1 &= b \\
\end{align}

引き算のしかたを具体的に書くと以下の通りです。 $x_k, y_k$ の漸化式は、 $q$ を拝借している点（および初期値）を除き $d_k$ と同じとなります。

\begin{alignat}{5}
& a \cdot x_{k-2} && + b \cdot y_{k-2} && = d_{k-2} \\
& a \cdot x_{k-1} && + b \cdot y_{k-1} && = d_{k-1} \\
&&& \Bigg \downarrow \quad q_{k-1} = \left \lfloor d_{k-2} / d_{k-1} \right \rfloor \\
& a \cdot \left( x_{k-2} - x_{k-1} q_{k-1} \right) && + b \cdot \left( y_{k-2} - y_{k-1} q_{k-1} \right) && = d_{k-2} - d_{k-1} q_{k-1} \\
\implies \quad & a \cdot x_k && + b \cdot y_k && = d_k \\
\end{alignat}

通常の互除法

（拡張する前の）ユークリッドの互除法は、2つの整数の最大公約数を効率よく求められるアルゴリズムです。

原理

「自然数 $a$ と $b$ の最大公約数は、（ $a$ と） $b$ と $|a-b|$ の最大公約数に等しい」ということを利用して、考える数を小さくしていくよう計算します。 $a = b$ ならそれが最大公約数です。

理屈がわかっていれば再帰関数で実装するのが手っ取り早いですが、手計算に倣ってループで実装することを目指します。

手計算

例として 2022 と 1224 の最大公約数を求めてみます。大きいほうから小さいほうを引くことを繰り返すと、

a		b
2022	＞	1224
798	＜	1224
798	＞	426
372	＜	426
372	＞	54
318	＞	54
...	＞	...
102	＞	54
48	＜	54
48	＞	6
42	＞	6
...	＞	...
12	＞	6
6	＝	6

…と長くなるので、普通は割り算を使います。例えば 372 から 54 を引き続けるのは、 $372 \div 54 = 6 余り 48$ なので、 6 回引くことによって 48 となります。

割り算を使った以下の表では、 $b$ と $r$ を次行の $a$ と $b$ にスライドさせて、 $a = b q + r$ という形を保っています。

a	b	商 q	余 r
2022	1224	1	798
1224	798	1	426
798	426	1	372
426	372	1	54
372	54	6	48
54	48	1	6
48	6	8	0
6	0	-	-

ゼロが出ると割り切れたことになり、直前の数値が最大公約数です。

実装

手計算での方法をコード化すると以下のようになります。大抵のプログラミング言語には剰余演算がありますが、今回は商と余りの関係も見れるように剰余演算なしで実装しています。

def gcd(a, b):
    while b != 0:
        q = a // b
        r = a - b * q  # == a % b
        a, b = b, r

    return abs(a)

※ 引数が負の整数の場合にも対応できるよう、最後に絶対値をとっています。

拡張互除法

解きたい問題

整数 $a, b$ が与えられたときに、次の方程式を満たす整数解 $x, y$ が存在します（ベズーの等式）。

a x + b y = \gcd (a, b)

※ なお、右辺を整数倍しても解があります（元の解を整数倍すればいいです）が、整数倍以外だと解が存在しません。したがって $a, b$ が互いに素（最大公約数が 1 ）のときに限り右辺がどんな整数でも解が存在します。

具体例として先ほどと同じ数をあてると以下のようになります。

2022 x + 1224 y = 6

しかし $x, y$ を実際に1組求めろと急に言われると、どこから手を付けていいのかさっぱりです。

手計算

右辺の条件は一旦忘れて、 $x, y$ に簡単な値を入れてみます。

\begin{align}
2022 \cdot 1 + 1224 \cdot 0 &= 2022 \\
2022 \cdot 0 + 1224 \cdot 1 &= 1224 \\
\end{align}

すると、2つの式を引き算したら、より小さな右辺を実現する $x, y$ を入れた式が出てきます。

\begin{alignat}{5}
& 2022 \cdot (1-0) && + 1224 \cdot (0-1) && = 2022-1224 \\
\implies \quad & 2022 \cdot 1 && + 1224 \cdot (-1) && = 798 \\
\end{alignat}

これを繰り返してみます。（もう一度最初の式から書いておきます）

\begin{alignat}{5}
& 2022 \cdot 1 && + 1224 \cdot 0 && = 2022 \\
& 2022 \cdot 0 && + 1224 \cdot 1 && = 1224 \\
& 2022 \cdot 1 && + 1224 \cdot (-1) && = 798 \\
& 2022 \cdot (-1) && + 1224 \cdot 2 && = 426 \\
& 2022 \cdot 2 && + 1224 \cdot (-3) && = 372 \\
& 2022 \cdot (-3) && + 1224 \cdot 5 && = 54 \\
\end{alignat}

右辺の変化に見覚えがあります。これは通常の互除法で最大公約数を求めたときと全く同じです。ならば引き算を続ければ、右辺が最大公約数になるときが訪れ、当初の方程式にたどり着くはずです。

この次を 372 - 54 - 54 - … と繰り返すのは面倒なので、何回引けばいいのかを割り算で求めてしまいます。通常の互除法でも扱った通り、商は 6 です。

\begin{alignat}{5}
& 2022 \cdot (2 - (-3) \cdot 6) && + 1224 \cdot ((-3) - 5 \cdot 6) && = 372 - 54 \cdot 6 \\
\implies \quad & 2022 \cdot 20 && + 1224 \cdot (-33) && = 48 \\
\end{alignat}

最後までやっていきます。

\begin{alignat}{5}
& 2022 \cdot (-3) && + 1224 \cdot 5 && = 54 \\
& 2022 \cdot 20 && + 1224 \cdot (-33) && = 48 \\
& 2022 \cdot (-23) && + 1224 \cdot 38 && = 6 \\
& 2022 \cdot 204 && + 1224 \cdot (-337) && = 0 \\
\end{alignat}

右辺がゼロになる直前が最大公約数であり、そのときの左辺にあるのが当初求めたかった解 $x, y$ のひとつです。また、最終行の式は他の式に何回加えても右辺を変化させないため、一般解は $(x, y) = (-23 + 204 k, 38 - 337 k)$ （ $k$ は任意の整数）と表せます。

実装

値の変化を表でも纏めておきます。

x	y	d	商 q	余 r
1	0	a = 2022
0	1	b = 1224	1	798
1	-1	798	1	426
-1	2	426	1	372
2	-3	372	1	54
-3	5	54	6	48
20	-33	48	1	6
-23	38	6	8	0
204	-337	0	-	-

d の列が通常の互除法で、その際に求まる q を拝借すれば同じ計算で x, y の列も求まります。

したがってコード化すると以下のようになります。

def ext_gcd(a, b):
    # a * x + b * y == d
    x1, y1, d1 = 1, 0, a
    x2, y2, d2 = 0, 1, b

    while d2 != 0:
        q = d1 // d2
        # r = d1 - d2 * q  # == d1 % d2
        d1, d2 = d2, d1 - d2 * q
        x1, x2 = x2, x1 - x2 * q
        y1, y2 = y2, y1 - y2 * q

    if d1 < 0:
        d1 = -d1
        x1, y1 = -x1, -y1

    return (d1, x1, y1)

※ 引数が負の整数の場合、求まる公約数 d1 が負になってしまう可能性があります。最大公約数には非負整数を選ぶべきなので、全体を $(-1)$ 倍して補正します。

応用：モジュラ逆数の計算

モジュラ逆数とは、ある整数 $a$ に対して合同式 $a x \equiv 1 \pmod M$ を満たす整数 $x$ のことです。分数の割り算が「逆数の掛け算」でできるように、合同式での割り算は「モジュラ逆数の掛け算」を考えれば楽に実行できるため重宝します。

合同式をやめて $a x + M y = 1$ と書き直すと、拡張互除法で解ける方程式に似ていることがわかります。大きな違いとして右辺が 1 に固定されているので、 $\gcd(a, M) = 1$ を満たさないと解が存在しません。 $a$ も $M$ も素数である必要はありませんが、法 $M$ が素数であれば任意の $a \not \equiv 0$ に対してモジュラ逆数 $x \equiv a^{-1}$ が求まることが保証されます。

def mod_inv(a, m):
    g, x, y = ext_gcd(a, m)

    if g != 1:
        raise Exception("a and m are not coprime.")

    if x < 0:
        x += m

    return x