0.はじめに

この記事の目的

ABC186-Eの中で用いる拡張ユークリッドの互除法の非再帰的実装が理解できなかったため、調べて理解した内容をこの記事で整理します。
拡張ユークリッドの互除法自体は、拡張ユークリッドの互除法〜一次不定方程式 ax + by = c の解き方〜など分かりやすい記事がたくさんあったため詳細は割愛し、この記事では、「1000000007 で割ったあまり」の求め方を総特集！〜逆元から離散対数まで〜に乗っているような拡張ユークリッドの互除法の非再帰実装の原理についてまとめます。

1.拡張ユークリッドの互除法の行列表現

いきなりですが、拡張ユークリッドの互除法を行列の形で表現してみます。

ax + by = \gcd(a,b)

を解く際には、

a = \lfloor\frac{a}{b}\rfloor b + r

として、$(a,b)$の問題を$(b,r)$の問題に置き換えて解いています。
これを行列として表現すると、

\begin{bmatrix}
a \\
b
\end{bmatrix}
=
\begin{bmatrix}
k & 1 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
b  \\
r
\end{bmatrix}

となります。（$\lfloor\frac{a}{b}\rfloor = k$ としています）。
拡張ユークリッドの互除法では、これを$(\gcd(a,b),0)$になるまで繰り返し行います。
そのため、$a=r_0,b=r_1,r=r_2, \dots,\gcd(a,b)=r_h$として、$k_i = \lfloor\frac{r_{i-1}}{r_i}\rfloor$と書くと、下記のように表現できます。

\begin{bmatrix}
r_0 \\
r_1
\end{bmatrix}
=
\begin{bmatrix}
k_0 & 1 \\
1 & 0
\end{bmatrix}
\dots
\begin{bmatrix}
k_{h-1} & 1 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
r_h  \\
0
\end{bmatrix} \tag{1}

ここで、

\begin{bmatrix}
k_i & 1 \\
1 & 0
\end{bmatrix}

の逆行列は、

\begin{bmatrix}
0 & 1 \\
1 & -k_i
\end{bmatrix}

なので、$(1)$式は、下記のように変形できます。

\begin{bmatrix}
r_h \\
0
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
1 & -k_{h-1}
\end{bmatrix}
\dots
\begin{bmatrix}
0 & 1 \\
1 & -k_0
\end{bmatrix}
\begin{bmatrix}
r_0  \\
r_1
\end{bmatrix} \tag{2}

最後に、

\begin{bmatrix}
x & y \\
u & v
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
1 & -k_{h-1}
\end{bmatrix}
\dots
\begin{bmatrix}
0 & 1 \\
1 & -k_0
\end{bmatrix} \tag{3}

と置けば、$(2)$式が

\begin{bmatrix}
r_h \\
0
\end{bmatrix}
=
\begin{bmatrix}
x & y \\
u & v
\end{bmatrix}
\begin{bmatrix}
r_0  \\
r_1
\end{bmatrix}

のように変形でき、この$(x,y)$が冒頭の$ax + by = \gcd(a,b)$の解$(x,y)$になります。

2.拡張ユークリッドの互除法の非再帰的表現

上の行列表現から、実際に$(x,y)$を求めるアルゴリズムを考えます。
$(3)$式をもとに計算すればいいので、初期値を$(x_0=1,y_0=0,u_0=0,v_0=1)$(単位行列)に設定します。

\begin{bmatrix}
x_1 & y_1 \\
u_1 & v_1
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
1 & -k_0
\end{bmatrix}
\begin{bmatrix}
x_0 & y_0 \\
u_0 & v_0
\end{bmatrix}

を計算すると、

\begin{array}{llll}
x_1 = u_0 \\
y_1 = v_0 \\
u_1 = x_0 - k_0u_0 \\
v_1 = y_0 - k_0v_0
\end{array}

になるので、プログラムとしては、

x -= k[0]*u
y -= k[0]*v
x,u = u,x
y,v = v,y

のようなイメージになります。
あとはこれを繰り返せばいいので、

def extGCD(a,b):
    x, y, u, v = 1, 0, 0, 1
    while b:
        k = a // b
        x -= k * u
        y -= k * v
        x, u = u, x
        y, v = v, y
        a, b = b, a % b
    return x, y

このように、非再帰の形で拡張ユークリッドの互除法を実装できます。

3.モジュラ逆数のみを求めたい場合(例.ABC186-E)

冒頭に書いたABC186-Eのように、競プロの場合は拡張ユークリッドの互除法はモジュラ逆数を求めるために使うことが多いと思います。
$(a,b)$が互いに素の場合に、$ax + by = 1$の解$x$が$ax \equiv 1(\mod b)$となるので、モジュラ逆数の算出に拡張ユークリッドの互除法が使えます。

この場合は、2.の最後のプログラムの中の$(v,y)$の値は特に必要ないので、下記のようなプログラムで算出することができます。
($x$を正の値にするために、ループの後の処理が一部追加されています。)

def modinv(a, b):
    p = b
    x, y, u, v = 1, 0, 0, 1
    while b:
        k = a // b
        x -= k * u
        y -= k * v
        x, u = u, x
        y, v = v, y
        a, b = b, a % b
    x %= p
    if x < 0:
        x += p
    return x