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

# はじめに
けんちょんさんの[「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜](https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a)という記事の[3-5. 拡張Euclidの互除法の逆元計算](https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a#3-5-%E6%8B%A1%E5%BC%B5-euclid-%E3%81%AE%E4%BA%92%E9%99%A4%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8B%E9%80%86%E5%85%83%E8%A8%88%E7%AE%97)において、modinv関数がなぜ以下のようになるのかを理解するのに時間がかかったので、解説を備忘録として残します。

https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a#3-5-%E6%8B%A1%E5%BC%B5-euclid-%E3%81%AE%E4%BA%92%E9%99%A4%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8B%E9%80%86%E5%85%83%E8%A8%88%E7%AE%97

```cpp: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$の逆元を求めます。
まず、以下の式を見てください
```math
ax + by = gcd(a, b)
```
今、$b$はmodinv関数中の入力`m`に対応し$1000000007$などの素数となるので、最大公約数は1となり
```math
ax + my = 1
```
となります。
次に、両辺の法を$m$とすると
```math
ax = 1 \mod m
```
となり、$m$を法とした$a$の逆元を求めることができます。

# modinv関数中の変数`u,v`は何を意味しているの？
自分が、理解するのに時間がかかった点として、modinv関数中の変数`u,v`の役割があります。この章では、これらがどうして必要でどんな役割を果たしているかを解説していきます。
## 準備
modinv関数は、拡張されたユークリッドの互除法に大きく寄ってます。
`u,v`の役割を知るために、まずは以下のwikipediaの「拡張された互除法」のセクションを読んでください。
[ユークリッドの互除法 - Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%A6%E3%83%BC%E3%82%AF%E3%83%AA%E3%83%83%E3%83%89%E3%81%AE%E4%BA%92%E9%99%A4%E6%B3%95#%E6%8B%A1%E5%BC%B5%E3%81%95%E3%82%8C%E3%81%9F%E4%BA%92%E9%99%A4%E6%B3%95)

## `u,v`が対応しているもの
準備より、等式 $mx+ny=gcd(m,n)$(今回の場合だと$ax+my=1$)の$x,y$において
```math
\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`に対応付けると

```math
\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$を掛けると
```math
\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行目における
```cpp
u -= t * v; swap(u, v);
```
を解説します。
これは、行列式で表すと
```math
\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関数は適切に逆元を求められていることが分かりました。

```cpp: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 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜](https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a)
・[ユークリッドの互除法 - Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%A6%E3%83%BC%E3%82%AF%E3%83%AA%E3%83%83%E3%83%89%E3%81%AE%E4%BA%92%E9%99%A4%E6%B3%95#%E6%8B%A1%E5%BC%B5%E3%81%95%E3%82%8C%E3%81%9F%E4%BA%92%E9%99%A4%E6%B3%95)