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

2020年3月21日2020年9月2日数学最大公約数,ユークリッドの互除法,拡張ユークリッドの互除法,1次方程式

最大公約数を求める高速なアルゴリズムとしてユークリッドの互除法が知られています。このユークリッドの互除法を拡張することにより、

$$ ax+by=gcd(a,b)$$

の形をした、2変数の一次方程式の整数解を求めることができます。

拡張ユークリッドの互除法のアルゴリズム

\(ax+by=gcd(a,b) \) を解くアルゴリズムは以下のようになります。

※ ユークリッドの互除法を逆順にさかのぼりながら、その時の a, b に対する x, y を計算していくイメージです。

※ x と y は別の再帰関数の呼び出し内でも値を共有しています。

※ extend_gcd の返り値は a と b の最大公約数 d になります。

C++ での実装例

long long extend_gcd(long long a, long long b, long long &x, long long &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    } else {
        long long d = extend_gcd(b, a % b, y, x);
        y -= a / b * x;
        return d;
    }
}

拡張ユークリッドの互除法の気持ちと正当性

アルゴリズムを理解するために、具体例を見てから、数学的考察をします。

具体例

48, 32 の最大公約数を例として考えてみます。これらの最大公約数は 16 になります。

$$ 48x + 32y = 16$$

それでは拡張ユークリッドの互除法を用いて求めてみましょう。

1. extend_gcd(48,32,x,y)を呼ぶ
   1. extend_gcd(32,16,y,x)を呼ぶ
      1. extend_gcd(16,0,x,y)を呼ぶ
         • b が 0 なので、xを1, yを0にして、16を返す
      2. d を 16 とする
      3. x(=1) から 32/16*y(=0) を引く
      4. d を返す
   2. d を 16 とする
   3. y(=0) から48/32*x(=1) を引く
   4. d を返す

これで終了しました。 返り値は通常の gcd と同様に最大公約数になり、今回は16です。x は 1, y は -1 となって終了して、実際に

$$ 48\times 1 + 32\times (-1) = 16$$

が成り立つので、これらは解になっています。

数学的考察

なぜ先程のアルゴリズムで求まるのかを考えてみましょう。

まず、以下の2式を考えます。

$$ ax + by = gcd(a,b) \tag{3}$$

$$ a = qb + r \tag{4}$$

(\(q\) は商、\(r\) は余りです）

(4)を(3)に代入して整理すると、

$$ b(qx+y) + rx = gcd(a,b)$$

のようになりました。この形は、(3)の式とよく似ているので、再帰的に求めることができます。

実際に、\(s:=qx+y\), \(t:=x\) と置き直すと、

$$bs + rt = gcd(a,b) = gcd(b,r) \tag{5}$$

のようになります。仮に \(s\), \(t\) が求まったとすると、\(x = t\), \(y = s – qx\) から容易に \(x\) ,\(y\) を求めることができます。

(3)の \(a\), \(b\) は、(5)の \(b\), \(r\) を見てもらえば分かるように、再帰的な呼出しをするにつれて小さくなっていきます。
最終的に \(a\), \(b\) は、通常のユークリッドの互除法と同様に、最大公約数と0になるので、

$$ gcd(a,b) \cdot s’ + 0\cdot t’ = gcd(a,b)$$

のようになります。故に、\(s’ = 1\), \(t’=0\)と定めてやれば、再帰をさかのぼって \(x\), \(y\) が決定できます。