swap関数を利用したユークリッドの互除法～割った余りの、そのまた余りの…～

ついこの前まで、ポインタについてシリーズ記事で解説していました。具体的には下のような感じ。

ざっくりと言うと、与えられた2つの引数を入れ替えるというswap関数を定義しようとしたけど、ポインタを使わないと狙い通りの動作になってくれないよ～って話でした。

大概のC言語の解説書だと「ポインタを使ってswap関数を当初の狙い通りの関数にできました。めでたしめでたし。」ってところで終わってるかと思います。ただ、せっかく関数を作ったら使ってみたい！というのが人情でしょう（僕だけかもしれませんが）。

しかし、入門者の方にいきなりswap関数を応用できるようなプログラムを考えろというのはすごく難しい事かと思います。そこで、僕なりに考えてswap関数を応用できるようなプログラムを作りました。

ということで、ポインタ解説は前回記事までで終わったのですが、今回はその補足というか、応用というかをやっていきます。

ユークリッドの互除法
具体的なアルゴリズム
C言語による実装
ソースコードの解説
1. EuclidAlgorithm関数について
  1. <ここから蛇足>
  2. </ここまで蛇足>
2. main関数について
まとめ

ユークリッドの互除法

先ほど言った「swap関数を応用できるようなプログラム」とは何かというと、ユークリッドの互除法を実行するプログラムです。

ユークリッドの互除法とは簡単に言うと、2つの自然数の最大公約数を求めるというアルゴリズムのことです。

どういうことかと言えば、例えば、12と21の最大公約数は3ですが、12と21というように2つの自然数だけを与えられたという条件の下で、3を導き出すというアルゴリズムがユークリッドの互除法ということです。

今の例みたいに12と21みたいにあまり大きくない自然数同士の最大公約数を求めるだけであれば、直観的に分かったり、直観的には分からなくても、2つの自然数それぞれの約数をすべて書きだしたりして、正しい答えを求めることができます。

しかし、与えられた2つの自然数が14599と25829だった場合など、数値が大きくなった場合はそのような方法では難しいと分かるかと思います。

しかし、このような最大公約数を求める方法を数学的に考えると、とても効率よく求められるアルゴリズムを導けるのです。

効率的に2つの自然数の最大公約数を求める方法がユークリッドの互除法と呼ばれるアルゴリズムになります。

今回の記事は数学的な話ではなく、あくまでもポインタ解説の延長線上の話という位置づけなので、数学的な裏付けについては解説しませんが、ご了承くださいませ～。

具体的なアルゴリズム

ユークリッドの互除法がどんなものか理解していただいた所で、具体的にどうすれば2つの自然数の最大公約数を効率的に求められるのかを解説していきます。

与えられた2つの自然数をaとbとしますが、2つの自然数には大小関係があるはずです。そこで、与えられた2つの自然数の内の、大きい方をa、小さい方をbとおくことにします（つまり、a>bが成り立つとします）。

そのとき、ユークリッドの互除法は次のようなアルゴリズムになります。具体例をアルゴリズムの後ろに付けているので、次のアルゴリズムが抽象的すぎて分かりづらいという方は参考にしてみてください。

aをbで割って、その商をq、余りをrとする。
（つまり、\( a \div b = qb + r \)となるようなqとrを求めます。分かりづらければ、下の解説を参照してください）
aをrに置き換える（つまり、2つの自然数の内の大きい方の数値を余りに更新する）。
bと新しいaを比較して、大きい方をaに、小さい方をbにする。
余りが0になるまで、ステップ1から3を繰り返す。
余りが0になったときのaが求めたい最大公約数になっている。

例えば、12と21が与えられた場合、この2つの自然数に対してユークリッドの互除法アルゴリズムを適用してみます。

まず、ユークリッドの互除法を適用するための準備として、2つの自然数の内の大きい方をa、小さい方をbとする必要があります。この場合、21がa、12がbになります。

そして、ここからがユークリッドの互除法の本体です。

ステップ1
aをbで割ったときの商qと余りrを求めます。この場合、\( 21=q \times 12 + r \)となるようなqとrを求めるので、q=1、r=9となります。q=0、r=21とかではないのでご注意を。あくまでもaをbで割るという割り算をしたときの商と余りなので、余りは0以上、12未満でないといけません。

ステップ2
r=9だったので、新しくa=r=9とします。

ステップ3
aとbを比較します。aはステップ2で9に変わっているので、今はa=9、b=12という状況です。なので、a<bとなってしまっています。aは2つの自然数の内の大きい方でなければいけないので、aとbを入れ替えて、a=12、b=9とします。

ステップ4
新しいaとb（つまり、12と9）に対して、余りが0（つまり、r=0）となるまで、ステップ1から3を繰り返します。当然、aとbもその都度更新していきます。具体的には、
12と9に対して → q=1、r=3（ステップ1）→ a=r=3（ステップ2）→ a=9、 b=3（ステップ3）
9と3に対して → q=3、r=0（ステップ1）→ 1=r=0（ステップ2）→ a=3、b=0（ステップ3）
という感じです。

ステップ5
ステップ4で、余りが0になるまでステップ1から3を繰り返した結果、最終的にはa=3、b=0という状況になっています。この3が求めたい最大公約数になります。

以上が、ユークリッドの互除法になります。

準備の段階とステップ3で、aとbの数値を入れ替えるという操作をしていましたが、この操作をするときにswap関数を使うことになります。

では、次の節で、ユークリッドの互除法をC言語で実装してみます。

C言語による実装

ソースコードは次の通りです。

/********************************
* ユークリッドの互除法により
* 最大公約数を求めるアルゴリズム
*********************************/

#include <stdio.h>
#include <limits.h>

void Swap( int *, int * );
int EuclidianAlgorithm( int, int );

int main()
{
  int x = 14599;
  int y = 25829;
  int gcd;

  printf( "x = %d, y = %d\n", x, y );

  gcd = EuclidianAlgorithm( x, y );

  printf( "%dと%dの最大公約数 = %d\n", x, y, gcd );

  return 0;
}

// 与えられた2つの引数を入れ替える
void Swap( int *i, int *j )
{
  int temp;

  temp = *i;
  *i = *j;
  *j = temp;

  return;
}

// ユークリッドの互除法により最大公約数を求める
// 求めた最大公約数を返す
int EuclidianAlgorithm( int a, int b )
{
  int surplus;

  printf( "ここから計算過程\n" );

  while( b != 0 )
  {
    // aの方が小さければ、aとbを入れ替える(でないと、後ろの処理が複雑になる)
    if( a < b ){ Swap( &a, &b ); }
    printf( "a = %d, b = %d\n", a, b );

    // 剰余を計算する
    // ただし、ゼロ除算を回避するためのif分岐が必要
    if ( ( b == 0 ) || ( ( a == LONG_MIN ) && ( b == -1 ) ) )
    { break; }
    else
    { surplus = a % b; }

    // ここではaの方がbよりも大きいことが保証されているので、
    // aをsuurplusに更新する
    a = surplus;
  }

  printf( "ここまで計算過程\n" );

  return a;
}

実行環境はMSYS2、gcc9.2.0です。こちらのプログラムを実行していただくと次のような結果になるかと思います。

x = 14599, y = 25829
ここから計算過程
a = 25829, b = 14599
a = 14599, b = 11230
a = 11230, b = 3369
a = 3369, b = 1123
a = 1123, b = 0
ここまで計算過程
14599と25829の最大公約数 = 1123

では、次の節でこのプログラムの解説を。

ソースコードの解説

swap関数についてはポインタ解説記事の中で説明していたので、ここではswap関数以外について説明します。swap関数についての説明が欲しい方は、ポインタ解説のシリーズ第1回目の記事をご覧くださいませ～（記事の中ほどにswap関数の動作が書いてあります）。

EuclidAlgorithm関数について

今回のメインテーマであるユークリッドの互除法はEuclidianAlgorithm関数で実行しています。この関数では2つの数値を引数として取っています。そして、その2つの数値の最大公約数を返り値としています。

EuclidianAlgorithm関数での中身ですが、whileループの中がユークリッドの互除法を実行する部分です。その中で説明の都合上、余りを格納しておくためだけの変数が欲しかったので、surplusというint型の変数を宣言しています。

まず、if文で引数として与えられた2つの数の大小関係を比較します。そして、もしaの方が大きければその後のユークリッドの互除法が使えないので、Swap関数でaとbを入れ替えています（準備、兼ステップ3）。

<ここから蛇足>

ここで注意していただきたいことがあります。それは、Swap関数が入れ替えた数値は何かということです。「え、それは引数として与えられたaとbじゃないの？」その通りです。その通りではあるのですが、それをメモリの話に置き換えるとどのようになるでしょうか？

もっと言えば、main関数のxとyは入れ替わるでしょうか？上のソースコードの実行結果を確認していただければ、EuclidianAlgorithm関数の後でxとyを呼び出してもその数値が変わっていないので、main関数のxとyは入れ替わらないと分かるのですが、それではSwap関数はどのような動作をしているのでしょうか？もし興味があれば考えてみてください。

</ここまで蛇足>

Swap関数で入れ替えるのは、if文が真のときのみです。言い換えれば、aがbよりも小さいときだけです。当然のことですが、aがbよりも大きいときは入れ替わらないのです。

つまり、このif文以降では、aやbに何かしらの操作が加わらない限り、aがbよりも大きい数値になっていることが保証されているのです。なので、この直後であれば、安心してユークリッドの互除法を実行できます。

プログラムのその次の処理はprintf関数ですが、これは単純に数値を表示しているだけです。ソースをさらに読み進めると、再びif文に行き当たります。具体的には次のようになっていました。

    if ( ( b == 0 ) || ( ( a == LONG_MIN ) && ( b == -1 ) ) )
    { break; }
    else
    { surplus = a % b; }

この部分は、パッと見はとても複雑なことをしてそうに見えますが、実際には剰余計算をしているだけです（ステップ1）。要するに、ただaをbで割った余りを求めているだけです。しかしそれなら、次のようにするだけでも良さそうです。

    surplus = a % b;

なぜ複雑なプログラムにしているのでしょうか？

簡単に言えば、エラー回避のためです。

C言語の規格ではX%Yと表記したとき、その第二引数（つまりこの場合はY）が0の場合、その動作は未定義とされています（C11の規格書の6.5.5 Multiplicative operatorsの第5項目に記載されています。ドラフトではありますが）。

そのため、bが0になってもらっていては困るわけです。他にもいくつかの注意点があるらしいのですが、それはこちらのページにまとめられています。今回使用している剰余演算のif文も同ページのものを使わせていただきました。

剰余演算が終わったら、余りをaに代入しています（ステップ2）。今回は分かりやすくするためにsurplusという剰余を保存しておくためだけの変数を用意しましたが、実際にはaに直に余りを代入した方がソースコードが短く済みますし、使用メモリも少なくできるので、余りはaに直に代入した方がいいと思います。

あとは、これをbが0になるまで繰り返すだけです（ステップ4）。繰り返すと、先に説明したSwap関数が実行されますが、これがステップ3になります。

そして、bが0になるとエラー回避用のif文でbreakが実行されます（ステップ4のループ抜け条件）。

ループを抜けた時点でaに入っている数値は最大公約数になっていますから、その数値をreturn文で返しています（ステップ5）。