Mister雑記

競プロやります。

「拡張ユークリッドの互除法」の実装

長いこと「拡張ユークリッドの互除法」の実装を飲み込めていないため、いい加減理解しようと思い立ってパパっと書きました。

内容はあくまで実装がメインなので、理論的な話はけんちょんさんの記事を読むのが一番確実かと思います。

また本記事は「自分にとっての分かりやすさ」を最優先に書いているので、普段より分かりにくい記事になっていると思われます。ご了承下さい。

目次

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

そもそも「拡張ユークリッドの互除法」とは、 a, b \in \mathbb{N}に対して、

  •  gcd(a, b)
  •  ax + by = gcd(a, b)なる x, y \in \mathbb{Z}

を両方とも求めるアルゴリズムを指す。ここで gcd(a, b)は「 a bの最大公約数(Greatest Common Divisor)」である。


1つ目は「ユークリッドの互除法」そのものである。求め方もユークリッドの互除法とほとんど同じとなる。

2つ目については、条件を満たす x, y1つだけ返す。これは

  • 解が1つ求まればそこから一般解を求めるのはそこまで大変ではない
  • 競プロでは1つ求まれば十分なケースが多い

といった都合上、機能としては十分である。

ユークリッドの互除法

「拡張」ユークリッドの互除法なので、このアルゴリズムを実装するにはまず「ユークリッドの互除法」を十分に理解する必要がある。

多くの場合、ユークリッドの互除法

  •  gcd(a, b) = gcd(b, a \bmod b)  (a \geq b)
  •  gcd(a, 0) = a

という性質を利用して、再帰関数で実装される。C++で実装すると以下の通り。

// aとbの最大公約数を返す
int gcd(int a, int b) {
    if (a < b) return gcd(b, a);  // a >= bに統一する
    
    if (b > 0) {
        return gcd(b, a % b);  // 1つ目のケース
    } else {
        return a;  // 2つ目のケース 終端条件
    }
}

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

それではメインの「拡張ユークリッドの互除法」の実装について考える。関数名は「 extgcd」(extended gcd)が一般的なので、ここでもその名前で呼ぶことにする。

関数 extgcdの返り値だが、これは gcd(a, b)だけであることが多い。では x, yはどうやって返すのかと言うと、「参照渡し」によって内部で書き換えられることで実現される。つまり型だけ見るとこうなる。

// gcd(a, b)を返す
// ax + by = gcd(a, b)となるようにx, yを書き換える

int extgcd(int a, int b, int& x, int& y);

これについては「こうすると実装しやすい」という程度に考えておけばいいと思われる。

遷移

引数の遷移は先の gcdと同じで、 (a, b) \rightarrow (b, a \bmod b)となる。つまり再帰呼出しによって

  •  bx' + (a \bmod b)y' = gcd(b, a \bmod b)

なる x', y'が先に求まる。ここから「 ax + by = gcd(a, b)」なる x, yを求めるにはどうすれば良いだろうか?1


まず右辺は gcd(a, b) = gcd(b, a \bmod b)なので共通している。よって左辺を上手く変形してやればいい。

そこで、「 a = (a / b) \cdot b + (a \bmod b)」という関係式に着目する。これは割り算の商と余りの関係から導かれる。この式を ax + by = gcd(a, b)の左辺に代入すると、

  •  ax + by = ( (a / b) \cdot b + (a \bmod b) ) x + by = b \cdot ((a / b) x + y) + (a \bmod b) x

となる。左辺だけまとめると以下の通り。

  •  b \cdot ((a / b) x + y) + (a \bmod b) x
  •  bx' + (a \bmod b)y'

 b, a \bmod bの係数を比較すると、

  •  x = y'
  •  y = x' - (a / b) x

という関係式が出てくる。これに沿って x, yを書き換えてやればいい。


これをそのまま実装しようとすると以下のようになる。

    if (b > 0) {
        int g = extgcd(b, a % b, x, y);
        int prex = x;
        x = y;
        y = prex - (a / b) * x;
    } ...

一方、変数を渡す順番を変えることでより短く実装できる。

    if (b > 0) {
        int g = extgcd(b, a % b, y, x);
        y -= (a / b) * x;
    } ...

実装がスマートな反面、初見だととても理解しにくい。一長一短なので、個人の好みで選べばいいと思われる。

終端条件

終端条件も gcdと同様に b = 0となる。

では ax + 0 \cdot y = aとなる xは何だろうか?...言わずもがな x = 1である。 yは何でもいいのだが、特別な理由がない限りは y = 0がいいだろう。

    if (b > 0) {
        ...
    } else {
        x = 1, y = 0;
        return a;
    }

実装例

後は a \lt bの場合にswapさせる処理を加えればよいが、これは簡単なので説明は割愛。

以上の実装をまとめると、以下のようになる。

// gcd(a, b)を返す
// ax + by = gcd(a, b)となるようにx, yを書き換える
int extgcd(int a, int b, int& x, int& y) {

    if (a < b) return extgcd(b, a, y, x);
    // xとyを入れ替え忘れないように注意
    
    if (b > 0) {
        int g = extgcd(b, a % b, y, x);
        y -= (a / b) * x;
        return g;
    } else {
        x = 1, y = 0;
        return a;
    }
}

最後に

拡張ユークリッドの互除法を実装できたのはいいのですが、実は今まで拡張ユークリッドの互除法を使う問題に出会ったことがありません。つまりverifyがまだできてないんですね。

ということでAtCoderCodeforces、yukicoderでverifyできる問題があったら教えていただけると幸いです。それと、もしもバグがあれば指摘していただけるととても助かります。


  1. ここから先の流れを自分で書けるようになれば、拡張ユークリッドの互除法は9割理解したと言っても過言ではない。そこまで飛躍した展開もないはずなので、根気よく追うことを推奨する。