逆元計算の解説 〜modinv関数を理解しよう〜

はじめに

けんちょんさんの「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜という記事の3-5. 拡張Euclidの互除法の逆元計算において、modinv関数がなぜ以下のようになるのかを理解するのに時間がかかったので、解説を備忘録として残します。

modinv関数

long long modinv(long long a, long long m) {
    long long b = m, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= m; 
    if (u < 0) u += m;
    return u;
}

考え方

ユークリッドの互除法を利用して$a$の逆元を求めます。
まず、以下の式を見てください

ax + by = gcd(a, b)

今、$b$はmodinv関数中の入力mに対応し$1000000007$などの素数となるので、最大公約数は1となり

ax + my = 1

となります。
次に、両辺の法を$m$とすると

ax = 1 \mod m

となり、$m$を法とした$a$の逆元を求めることができます。

modinv関数中の変数u,vは何を意味しているの？

自分が、理解するのに時間がかかった点として、modinv関数中の変数u,vの役割があります。この章では、これらがどうして必要でどんな役割を果たしているかを解説していきます。

準備

modinv関数は、拡張されたユークリッドの互除法に大きく寄ってます。
u,vの役割を知るために、まずは以下のwikipediaの「拡張された互除法」のセクションを読んでください。
ユークリッドの互除法 - Wikipedia

u,vが対応しているもの

準備より、等式 $mx+ny=gcd(m,n)$(今回の場合だと$ax+my=1$)の$x,y$において

\begin{pmatrix}
x & y \\
u & v
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
1 & -k_{h-1}
\end{pmatrix}
\begin{pmatrix}
0 & 1 \\
1 & -k_{h-2}
\end{pmatrix}
\cdots
\begin{pmatrix}
0 & 1 \\
1 & -k_{0}
\end{pmatrix}

が成り立つことが分かりました。
ここで、肝心のu,vが何なのかですが、先に答えを言うと、これらはそれぞれ上式の$x,u$に対応してます。
実際に、modinv関数中の変数u,vそしてtに対応付けると

\begin{pmatrix}
u & \_ \\
v & \_
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
1 & -t_{h-1}
\end{pmatrix}
\begin{pmatrix}
0 & 1 \\
1 & -t_{h-2}
\end{pmatrix}
\cdots
\begin{pmatrix}
0 & 1 \\
1 & -t_{0}
\end{pmatrix}

のように書き直せます。
さらにアンダーバーを消すために、両辺に右から$(1 \ 0)^T$を掛けると

\begin{pmatrix}
u \\
v
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
1 & -t_{h-1}
\end{pmatrix}
\begin{pmatrix}
0 & 1 \\
1 & -t_{h-2}
\end{pmatrix}
\cdots
\begin{pmatrix}
0 & 1 \\
1 & -t_{0}
\end{pmatrix}
\begin{pmatrix}
1 \\
0
\end{pmatrix} \quad (A)

となります。ここで$1$ $0$はそれぞれu,vの初期値に対応します。
以上より、(A)式をコードで再現できれば、$u$として$ax + my = 1$における$a$の逆元$x$が求められることが分かりました。

6行目u -= t * v; swap(u, v);の意味

最後にmodinv関数の6行目における

u -= t * v; swap(u, v);

を解説します。
これは、行列式で表すと

\begin{pmatrix}
u \\
v
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
1 & -t \\
0 & 1
\end{pmatrix}
\begin{pmatrix}
u \\
v
\end{pmatrix}　\\

\because \
\begin{pmatrix}
u \\
v
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
1 & -t
\end{pmatrix}
\begin{pmatrix}
u \\
v
\end{pmatrix}

となります。
もう、お気づきですね。そうです、u,vの初期値を$1,0$として、この処理を繰り返していくことで、最終的に(A)式になります。
よって、modinv関数は適切に逆元を求められていることが分かりました。

modinv関数再掲

long long modinv(long long a, long long m) {
    long long b = m, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= m; 
    if (u < 0) u += m;
    return u;
}

最後に

今回は、modinv関数の解説をしました。コードの理解に困っている人の助けになれば幸いです。
間違いなどあれば、気軽にご指摘ください。

参考

・「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜
・ユークリッドの互除法 - Wikipedia