超わかりやすい拡張ユークリッドの互除法

数a,bが与えられる。
1. 最初はx=a, y=bとする。
2. xをyで割った余りにする。
3. yをxで割った余りにする。
4. 2に戻る。
5. xかyが0になったときに、もう一方が最大公約数である。

ここで、この数x,yを、別の方法で表現することを考える。つまり

x = a * x' + b * x''
y = a * y' + b * y''

まず、最初の初期化では次のようになる。

x = a * 1 + b * 0 = a
y = a * 0 + b * 1 = b

また、この形式を保持しながら、「余り」の計算を行うことが可能である。
xをyで割った余りをmod、割った整数の商をdivとすると、

mod = x - y * div

が成立することを利用する。

         y = a *    y'    + b *    y''
=> div * y = a * (div*y') + b * (div*y'')

         x = a *    x'    + b *    x''
+) div * y = a * (div*y') + b * (div*y'')
-------------------------------------------------
       mod = a * (x'-div*y') + b * (x''-div*y'')

この計算を、a%bの計算のついでに行うだけで、
最終的には

gcm = a * something + b * something

という形の式が得られる。

というわけで、簡単な拡張ユークリッドのデモを書いた。

#include <cstdio>
#include <utility>
#include <algorithm>

using namespace std;

pair<int, pair<int, int> > ee(int a, int b) {
    pair<int, pair<int, int> > x(a,make_pair(1,0));
    pair<int, pair<int, int> > y(b,make_pair(0,1));
    int div;
    while(true) {
        if(!y.first)return x;
        div = x.first/y.first;
        x.first         = x.first         - div * y.first;
        x.second.first  = x.second.first  - div * y.second.first;
        x.second.second = x.second.second - div * y.second.second;
        if(!x.first)return y;
        div = y.first/x.first;
        y.first         = y.first         - div * x.first;
        y.second.first  = y.second.first  - div * x.second.first;
        y.second.second = y.second.second - div * x.second.second;
    }
}
// mod and num must be co-prime.
int inverse(int mod, int num) {
    pair<int, pair<int, int> > eret = ee(num, mod);
    return ((eret.second.first%mod)+mod)%mod;
}
// mod must be prime number.
int inverse2(int mod, int num) {
    int pw = mod-2;
    int ret = 1;
    for(int i = 30; i >= 0; i--) {
        ret = (ret*ret)%mod;
        if((pw>>i)&1) {
            ret = (ret*num)%mod;
        }
    }
    return ret;
}

int main(int argc, char **argv) {
    pair<int, pair<int, int> > eret = ee(936, 102);
    printf("936 * %d + 102 * %d = %d\n", eret.second.first, eret.second.second, eret.first);
    for(int i = 1; i < 47; i++) {
        int inv = inverse(47, i);
        int inv2 = inverse2(47, i);
        printf("inverse:  %2d * %2d === 1 (mod 47)\n", i, inv);
        printf("inverse2: %2d * %2d === 1 (mod 47)\n", i, inv2);
    }
    return 0;
}

剰余環における逆元を拡張ユークリッドで求めるのだが、同時にフェルマーの小定理とバイナリ法を利用した方法(こちらは素数を法とした場合しか計算できない)も示した。

他には、拡張ユークリッドを使うと、中国人剰余定理に応用できるので、これをさらに百五減算とか多倍長とかに応用できると思う。