拡張Euclid互除法による逆元の計算 の コード解説

はじめに

キッカケは「けんちょんさん」のMOD計算をするためのライブラリ modint の割り算を見た時に、「はぁわけわからん」と思ったこと

modint の全容とかは調べてもらうとして、対象のコードはこちら。modint ではなく、単純に逆元を求める方のコードを採用（modint の方も求めた逆元をかけて割り算するだけなので）。ちなみにリンクはこちら (記事内容：「「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜」「 3.5 拡張 Euclid の互除法による逆元計算」)

拡張ユークリッドの互除法による逆元を求めるコード１(modinv)

# include <iostream>
using namespace std;

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;
}

int main() {
    // mod. 13 での逆元を求めてみる
    for (int i = 1; i < 13; ++i) {
        cout << i << " 's inv: " << modinv(i, 13) << endl;
    }
}

この modinv が激しく謎だった。
全部読めば、このコードの意味がわかるよ。

前提

• 逆元は既に知っている
• ユークリッドの互除法も知っている
• 拡張ユークリッド互除法 ← はぁ？知らんがな

拡張ユークリッドの互除法による逆元の高速計算方法自体に関しては、けんちょんさんのこちらの記事を参照

けんちょんさんは神。
どうでもいいけど、「けんちょん」さんって呼び方は「さかなくん」さんと同じ感じがする

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

拡張ユークリッドの互除法 := 一次不定方程式 $ax + by = c$ の整数解 x, y を求めるアルゴリズム

ちなみに「不定方程式」とは、方程式の数よりも未知変数の数が多い方程式のこと

ではなぜ、拡張ユークリッドの互除法によって逆元を求められるかは 「「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜」「 3.5 拡張 Euclid の互除法による逆元計算」 のコードの前まで読むことでわかる。

実際の拡張ユークリッドの互除法による解の求め方は、けんちょんさんの記事 「拡張ユークリッドの互除法 〜 一次不定方程式 ax + by = c の解き方 〜」「3. 一次不定方程式 ax + by = c の整数解」 を参照しましょう。

整理すると、「拡張ユークリッドの互除法を使って、逆元を求める」→「mod p における a の逆元 x を求めるために、拡張ユークリッドの互除法を使う」
↓
$ ax + py = 1 $ の解(x, y)を拡張ユークリッドの互除法により解(x, y)を求める。
その解 x が a の逆元ということがわかる。pはmod。

これを素直にコードに落とし込むと、こちらのようになる。

以下、コード

拡張ユークリッドの互除法による逆元を求めるコード(extGCD)

// 返り値: a と b の最大公約数
// ax + by = gcd(a, b) を満たす (x, y) が格納される
long long extGCD(long long a, long long b, long long &x, long long &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    long long d = extGCD(b, a%b, y, x);
    y -= a/b * x;
    return d;
}

ここまでで「拡張ユークリッドの互除法を使って、逆元を求める」コードはわかった。
しかし、このコード自体は最初に紹介したものとは違う。
ようやく本題へ。

本題 (modinv関数がなぜうまくいくのかわからない)

まずは $111x + 30y = 3$ を例に拡張ユークリッドを追っていきましょう。

実際の解き方は先ほども紹介したこちら(「拡張ユークリッドの互除法 〜 一次不定方程式 ax + by = c の解き方 〜」「3-4.拡張ユークリッドの互除法のアルゴリズム的記述」)を参考。その他、変数なども準拠。
ちなみにこの記事における、ユークリッドの互除法の終了条件とは「$a = p * b + r$ において、$r = 0$ の時の b がgcd(a, b) = d」

まずは extGCD関数
以下の見方

1. 式を上から下に見ていって式変形を確認
   $a = p * b + r$ における b と r が次のa と b になる

2. 一番底から再帰的に求められる s, t を確認
   $x = t, y = s - qt$
   ここにおける、x, y とは1つ上の s, t のこと。繰り返すことで当初の x, y に行き着く

| 式 | 商 p || s | t | y = s - qt (s, t は一段下。y は t) |
|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
| $111x + 30y = 3$ |3| | 3 | -11 | = -2 - 3 * 3 |
| $30x + 21y = 3$ |1| | -2 | 3 | = 1 - 1 * (-2) |
| $21x + 9y = 3$ |2| | 1 | -2 | = 0 - 2 * 1 |
| $9x + 3y = 3$ |3| | 0 | 1 | = 1 - 3 * 0 |
| $3x + 0y = 3$ |-| | 1 | 0 | 初期値 |

まあ関数通りというかなるほど納得。

では次に、modinv関数について。 改めてコードを見ると、以下のことがわかる

拡張ユークリッドの互除法による逆元を求めるコード１(modinv)

# include <iostream>
using namespace std;

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;
}

long long b = m, u = 1, v = 0;

b が mod　つまり$ax + py = 1$で、(u, v)が(s, t)ってことだな

while (b) {

b が 0 になるまで、つまりユークリッドの互除法で $a = qb + r$ の b が 0 になるまで

long long t = a / b;

この t は商 p のこと

a -= t * b; swap(a, b);

$a = a - bp$ なので、つまり余りrのこと。それをswapしているので、(b, r)を次の(a, b)にしているということ。

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

$u = u - q * v$ つまり $y = s - qt$ それをswap してるので、extGCD関数でやったことと同じだな...**違うな？？？**なんで同時進行してんの？

再帰じゃないから、q の値が違ってうまくいくわけないじゃん！
疑問（本題）：なんでこのコードうまくいくの？？？

ここまで長かった。というわけで先ほど同様の図を使う。
今回はs, tが上から同時に求められることに注意。

| 式 | 商 p || s | t | y = s - qt (s, t, q は一段上。y は t) |
|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
| $111x + 30y = 3$ |3| | 1 | 0 | 初期値 |
| $30x + 21y = 3$ |1| | 0 | 1 | = 1 - 0 * 3 |
| $21x + 9y = 3$ |2| | 1 | -1 | = 0 - 1 * 1 |
| $9x + 3y = 3$ |3| | -1 | 3 | = 1 - (-1) * 2 |
| $3x + 0y = 3$ |-| | 3 | -10 | = -1 - (3) * 3 |

求めたい答え（s の最終値 3)は一致する。t は一致する。

ここで、$x = t$なので、extGCD関数では下から上に(s, t)の値を求めたので tの値が左上に、modint関数では上から下に求めたので tの値が左下に行っていることを感じる。
そうすると、s の値は、必ず t の値から来ているので、漸化式っぽいなと思う。

ここで改めて、extGCD に関して、それぞれの変数に式番号として、番号を振ってみる

| 番号 | 式 | 商 p || s | t | y = s - qt (s, t は一段下。y は t) |
|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
| 4 | $111x + 30y = 3$ |3| | 3 | -11 | = -2 - 3 * 3 |
| 3 | $30x + 21y = 3$ |1| | -2 | 3 | = 1 - 1 * (-2) |
| 2 | $21x + 9y = 3$ |2| | 1 | -2 | = 0 - 2 * 1 |
| 1 | $9x + 3y = 3$ |3| | 0 | 1 | = 1 - 3 * 0 |
| 0 | $3x + 0y = 3$ |-| | 1 | 0 | 初期値 |

そうすると $t_0 = t;(= 0),;; t_1 = s;(= 1)$ と見える。
つまり、一般式は $t_i = t_{i-2} - q_i * t_{i-1}$
なので、$t_2 = t_0 - q_2 * t_1$, $t_3 = t_1 - q_3 * t_2$ となる。

今回の s は $t_3$ であるので、$t_2$の式を代入して、整理すると以下のようになる。

$t_3 = (1 + q_3q_2)s - q_3t$

つまり、t = 0 なので、最終的な解 x を求めるのに必要なのは、一番最初と最後以外の商の積。
これは順不同なので modint 関数のようは 再帰をしなくてもその時々の商を掛け算しとけば上手くいくってことなんでしょう。多分、知らんけど。

二項係数との関連について

二項係数 nCk を求める時、逆元が必要なのは周知の事実。例の如く

では 逆元(inv)を求める時にどうするか
よくやる二項係数 (nCk mod. p)、逆元 (a^-1 mod. p) の求め方が大変ためになる。

今回やったのは、「1-3. 拡張 Euclid の互除法による逆元計算」
上記の記事で見た方がいいのは、「1-4. 改めて、inv[1], inv[2], ..., inv[n] を高速に計算する方法」「1-5. 逆元漸化式のもう 1 つの導出方法」

おわりに

まぁライブラリだし！！中身知らんくても使えればいいよ！！　納得できないものはあまり使いたくない質だから考えただけ。
間違ってても知らん。