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

数学

2つの整数の最大公約数を求める際、単純に素因数分解して、共通部分を求めようとすると時間がかかってしまいます。

このような時に、対数のオーダーで高速に最大公約数を求めるアルゴリズムとして、ユークリッドの互除法が古くから知られています。

さらに、このユークリッドの互除法を拡張することにより、

$$ax+by=c$$

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

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

まずはアルゴリズムを説明して、具体例を見てから、数学的証明を与えます。

原理としては単純で、「 a と b の最大公約数」と「 b と b割るaの余り の最大公約数」が等しくなることを利用しています。

これを利用して、a と b の最大公約数をもとめる gcd(a,b) のアルゴリズムは以下のようになります。

  1. b が 0 ならば、a を返す
  2. b が 0 でなければ、gcd(b, b%a) を返す

a,bの大小関係は考慮していません。最初の入力が b>a であったとしても、一度gcdを再帰的に呼び出したら、gcd(b,a)と順番が入れ替わることに注意してください。

具体例

48, 32 の最大公約数を例として考えてみます。まずは素因数分解を用いて求めてみると、

$$ 48 = 2^4\times 3$$

$$ 32 = 2 ^5$$

となるので、最大公約数は 16 になります。

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

  1. gcd(48,32) = gcd(32, 48%32)
  2. gcd(32,16) = gcd(16, 32%16)
  3. gcd(16,0) = 16

これで 16 が求まりました。

計算量

ユークリッドの互除法をのステップを2回行ったあとを考えると、gcd(a,b)のa,bはもともとの半分以下になっているはずです。

よって2回ごとに最悪でも半分になることを考えると、計算量は\(O(log ~min(a,b))\)になります。

数学的証明

以下の式が成り立つ事を示しましょう。

$$gcd(a,b) = gcd(b,b\%a)$$

これは、

$$gcd(a,b) \le gcd(b,b\%a)\tag{1}$$

$$ gcd(a,b) \ge gcd(b,b\%a) \tag{2}$$

の2つを証明することで示すことができます。

以下、\(a \ge b\)として考え、

$$ a = qb + r $$

のように表すことができるとします。(qは商、rは余りです)。

まず(1)を考えます。
a,b の任意の公約数をmとします。これが b, b%a の公約数でもあれば 、(1)の式が成立するはずです。

実際に、b は m の倍数で、

$$ b\%a = r = qb – a$$

から b%a も m の倍数です。

次に(2)を考えます。
(1)の逆で、b, b%a の任意の公約数をmとします。これが a, b の公約数でもあれば 、(2)の式が成立するはずです。

実際に、b は m の倍数で、

$$ a = qb +r $$

から、a も m の倍数です。

以上より示されました。

プログラム例

C++の例

long long gcd(long long a, long long b) {
    if (b == 0) {
        return a;
    } else {
        return gcd(b, a % b);
    }
}

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

先程のユークリッドの互除法のアルゴリズムを拡張することで、

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

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

最大公約数を返し、x,y の値を決める extend_gcd(a,b,x,y) は以下のアルゴリズムです。

  1. b が 0 ならば、x を 1、y を 0 にして、a を返す
  2. b が 0 でなければ、
    1. extend_gcd(b, a%b, y, x) を再帰的に呼び出して結果をdとする
    2. y から a/b*x を引く
    3. 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, t を見てもらえば分かるように、再帰的な呼出しをするにつれて小さくなっていきます。
最終的に a,b は、通常のユークリッドの互除法と同様に、最大公約数と0になるので、

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

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

プログラム例

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

数学

Posted by cs-learner