ユークリッドの互除法

最大公約数と最小公倍数

最大公約数（ＧＣＤ：Greatest Common Divisor）とは、２つの正整数ａとｂを素因数分解してときの共通項です。例えば
　　　ａ＝９２４（＝２×２×３×７×１１）
　　　ｂ＝３６０（＝２×２×２×３×３×５）
の最大公約数Ｇは
　　　Ｇ＝２×２×３＝１２
になります。

なお、最小公倍数（ＬＣＭ：Least Common Multiple）は、ａとｂで割り切れる最小の整数のことです。
　　　ａ＝Ｇ×（７×１１）
　　　ｂ＝Ｇ×（２×３×５）
なので、最小公倍数Ｌは、
　　　Ｌ＝Ｇ×（７×１１）×（２×３×５）
となります。すなわち、
　　　Ｌ＝ａ×ｂ／Ｇ＝９２４×３６０／１２＝２７７２０
で計算できます。

ユークリッドの互除法

ユークリッドの互除法は、最大公約数を計算する効率的なアルゴリズムです。

• アルゴリズム
  正整数ａ、ｂ（ａ＞ｂ）の最大公約数を求める手順を示します。
  　　ア　ａをｂで割った剰余をｒとする。
  　　イ　ｒが０ならばｂが最大公約数である。
  　　ウ　ｒが０でないならば、ａ←ｂ、ｒ←ｂとして、アへ戻る。
• プログラム
  function euclid(a, b) {
  　　r = a % b;　　　　/* ｒはａ／ｂの剰余 Ｃ言語なら a mod b と記述 */
  　　while (r > 0) {
  　　　　a = b;
  　　　　b = r;
  　　　　r = a % b;
  　　}
  　　return b;
  }
• 数値例（ａ＝９２４、ｂ＝３６０）
  　ａ　　　ｂ　　　ｒ
  ９２４　３６０　２０４（＝９２４－２×３６０）　　　（ア）
  ３６０　２０４　１５６（＝３６０－１×２０４）
  ２０４　１５６　　４８（＝２０４－１×１５６）
  １５６　　４８　　１２（＝１５６－３×４８）
  　４８　　１２　　　０（＝４８－４×１２）
  　　　　　 └─→最小公倍数＝１２

再帰関数によるユークリッドの互除法

ユークリッドの互除法は、再帰関数を用いると簡潔な表現になります。

• アルゴリズム
  　　　　　　　　　┌ a　　　　　　　　　　（b＝0 のとき）
  　　　gcd(a,b) ＝ ┤
  　　　　　　　　　└ gcd(b, a mod b)　　　（b＞0 のとき）
• プログラム
  　　function gcd(a, b) {
  　　　　if (b > 0) return( gcd(b, a % b) );
  　　　　else return(a);
  　　}
• 数値例（ａ＝９２４、ｂ＝３６０）
  　　　　 b　　a　　 b　　 a % b→b
  　　gcd(360, 924 % 360)　　204
  　　gcd(204, 360 % 204)　　156
  　　gcd(156, 204 % 156)　　 48
  　　gcd( 48, 156 % 48)　　 12
  　　gcd( 12, 48 % 12)　　　0
  　　　　└─→最小公倍数＝１２
　　　　

（注）プログラムの改善

上のプログラムには、いくつかの改善点があります。

• 引数として渡されたａとｂの値が変化してしまいます。呼び出したプログラムのなかでａやｂを用いることもあるので、一般的には、呼び出された側で新しい変数にして、影響を少なくするようにします。ここでは、引き渡される引数の名称を変えました。
• ａ＞ｂ＞０の制約がありますが、渡された数値の大きいほうをａとするように内部で置換したほうが、取り扱いを容易になります。また、ａやｂが負数のときにも使えるように、絶対値にしたほうがよいでしょう。
• 剰余ｒを求めた結果が、ｒ＞ｂ－ｒ　のときは、ｒ＝ｂ－ｒと置き換えたほうが、ｒの値が小さくなるので、より容易に計算することができます。
  　例えば、上のアでは、ｂ＝３６０、ｒ＝２０４ですので、ｒ＝３６０－２０４＝４８とするのです。
  　これは、９２４－３×３６０＝－４８からも求められます。後述の「証明」からわかるように、整数商ｑは最大公約数に影響を与えないので、このようにしてもよいのです。

これらを考慮したプログラムを示します。
　　function euclid(n1, n2) {
　　　　var a, b, r;　　/* これらの変数は関数内でのみ有効。他のプログラムには無影響 */
　　　　a = Math.abs(n1);　　/* n1の絶対値 */
　　　　b = Math.abs(n2);
　　　　if (b > a) {　　　　/* ｂ＞ａ のときの置き換え */
　　　　　　r = a;　a = b;　b = r;
　　　　}
　　　　r = a % b;
　　　　while (r > 0) {
　　　　　　if (r > b-r) r = b-r;　　/* 剰余の修正 */
　　　　　　a = b;
　　　　　　b = r;
　　　　　　r = a % b;
　　　　}
　　　　return b;
　　}

ユークリッドの互除法の証明

ａをｂで割ったときの商をｑ、剰余をｒとすると、
　　　ｒ＝ａ－ｑ×ｂ
となります（ａ＝９２４、ｂ＝３６０とすると、ｒ＝９２４－２×３６０＝２０４）。
ここで、ａとｂの最大公約数 ＧＣＤ(ａ,ｂ) をｇとすれば、ａもｂもｇで割り切れるので、ｒもｇで割り切れます（この段階では、ｇが２４であることはわかっていないのですが、ｒ＝２０４（＝１７×２４）はｇ＝２４で割り切れます）。すなわち、
　　　ＧＣＤ(ａ,ｂ)＝ＧＣＤ(ｂ,ｒ)　　　ＧＣＤ(９２４,３６０)＝ＧＣＤ(３６０,２０４)
の関係があります。ＧＣＤ(ａ,ｂ) を求めるかわりに、ＧＣＤ(ｂ,ｒ) を求めればよいし、しかも、ｒ＜ｂ＜ａなので、ＧＣＤ(ｂ,ｒ) を求めるほうが簡単です。
「ｂ←ａ、ｒ←ａ／ｂの剰余」の置き換えを繰り返すことにより、より小さい数の組み合わせになります。 そして、ｒが０になるまで繰り返されるのですから、最終的には、
　　　ＧＣＤ(ａ,ｂ)＝ＧＣＤ(ｂ,０)
となります（ＧＣＤ(３６０,２０４)＝ＧＣＤ(２０４,１５６)＝ＧＣＤ(１５６,４８)＝ＧＣＤ(４８,１２)＝ＧＣＤ(１２,０) ）。
ここで、ＧＣＤ(ｂ,０) とは、ｂと０との最大公約数であり、０はどんな数でも割り切れるのですから、ＧＣＤ(ｂ,０) ＝ｂ であることは明白です。すなわち、
　　　ＧＣＤ(ａ,ｂ)＝ｂ　　　ＧＣＤ(９２４,３６０)＝１２
となります。

計算量

最大公約数の計算量では、剰余の計算回数が尺度になります。私は理解していないのですが、ラメ（Lame'）の定理によれば、最悪の場合でも、剰余の計算回数ｂの桁数の５倍以下であることが証明されています（ｂ＝３６０は３桁、３×５＝１５回以下）。
　また、ｋ回の剰余計算が必要となるとき、ｂはｋ番目のフィボナッチ数以上であることが証明されています（フィボナッチ数とは、Ｆ(ｋ)＝Ｆ(ｋ－１)＋Ｆ(ｋ－２)、　Ｆ(１)＝１、Ｆ(２)＝１で定義される級数のこと。Ｆ(１３)＝２３３，Ｆ(１４)＝３７７）。

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

拡張ユークリッドの互除法は、自然数ａとｂを与えたとき、不定方程式
　　　ａｘ＋ｂｙ＝ｇ　（ｇはａとｂの最大公約数）
となる整数（負の場合もある）ｘとｙを求める方法です。

例えば、ａ＝９２４，ｂ＝３６０とすると、ｇ＝１２となり、
　　　９２４×（－７＋３０ｔ）＋３６０×（１８－７７ｔ）＝１２
になります。このような、
　　　ｘ＝－７＋３０ｔ
　　　ｙ＝１８－７７ｔ　　（ｔは任意の整数）
を求める方法です。

アルゴリズム

ａｘ＋ｂｙ＝ｇを
　　　ａ×ｘ_ｉ＋ｂ×ｙ_i＝ｒ_ｉ　　　①
と表記します。
　初期値を
　　　ｘ_０＝１，ｙ_０＝０、ｒ_０＝ａ（＝９２４）
　　　ｘ_１＝０，ｙ_１＝１、ｒ_１＝ｂ（＝３６０）
とすれば、①は
　ｉ＝０のとき、ａ×ｘ_０＋ｂ×ｙ_０＝ｒ_０　→　ａ×１＋ｂ×０＝ａ
　ｉ＝１のとき、ａ×ｘ_１＋ｂ×ｙ_１＝ｒ_１　→　ａ×０＋ｂ×１＝ｂ
となり、成立しています。

ｑ_ｉをｒ_ｉ－１／ｒ_ｉの整数部分、ｒ_ｉ＋１をその剰余とし、
　　　ｒ_ｉ＋１＝ｒ_ｉ－１－ｑ_ｉ×ｒ_ｉ　　②
　　　ｘ_ｉ＋１＝ｘ_ｉ－１－ｑ_ｉ×ｘ_ｉ　　③
　　　ｙ_ｉ＋１＝ｙ_ｉ－１－ｑ_ｉ×ｙ_ｉ　　④
と定義します。

ｉ＝１のときは、ｑ_１＝ｒ_０／ｒ_１の整数部分なので、ｑ_１＝９２４／３６０＝２となります。
②③④から、
　　　ｒ_２＝ｒ_０－ｑ_１×ｒ_１　（９２４－２×３６０＝２０４）
　　　ｘ_２＝ｘ_０－ｑ_１×ｘ_１　（１－２×０＝１）
　　　ｙ_２＝ｙ_０－ｑ_１×ｙ_１　（０－２×１＝－２）
となるので、これを①に代入すれば、
　　　ａ×ｘ_２＋ｂ×ｙ_２＝９２４×１－３６０×（－２）＝２０４＝ｒ_２
となり、①が成立しています。

ｉ＝ｋのときは
　　　ｘ_ｋ＋１＝ｘ_ｋ－１－ｑ_ｋ×ｘ_ｋ
　　　ｙ_ｋ＋１＝ｙ_ｋ－１－ｑ_ｋ×ｙ_ｋ
とすれば、
　　　ａ×ｘ_ｋ＋１＋ｂ×ｙ_ｋ＋１＝（ａ×ｘ_ｋ－１＋ｂ×ｙ_ｋ－１）＋ｑ_ｋ×（ａ×ｘ_ｋ＋ｂ×ｙ_ｋ）
となり、
　　　ａ×ｘ_ｋ－１＋ｂ×ｙ_ｋ－１＝ｒ_ｋ－１
　　　ａ×ｘ_ｋ＋ｂ×ｙ_ｋ＝ｒ_ｋ
ですから、
　　　ａ×ｘ_ｋ＋１＋ｂ×ｙ_ｋ＋１＝ｒ_ｋ＋１
となり、数学的帰納法により、①が常に成立することが証明されます。

ここで、ｒ_ｉ－１＝０になるときのｒ_ｉがａとｂの最大公約数であることは、先のユークリッドの互除法での説明からわかります。

そして、①の関係を保ちながら式の変形をしてきたので、ｘ_ｉとｙ_ｉは、
　　　ａｘ＋ｂｙ＝ｇ
を満足しています。
　後述のトレースにより、ａ＝９２４、ｂ＝３６０のときは、ｘ＝－７，ｙ＝１８になり、
　　　９２４×（－７）＋３６０×（１８）＝１２
となります。

不定方程式では、ｘとｙの値はいくらでもあります。証明は省略しますが、上の計算で得られたｘとｙをｘ^*、ｙをｙ^*とし、最大公約数をｇとしたとき、ｔを任意の整数として
　　　ｘ＝ｘ^*＋(ｂ／ｇ)×ｔ　　（ｘ＝－７＋３０ｔ）
　　　ｙ＝ｙ^*－(ａ／ｇ)×ｔ　　（ｘ＝１８－７７ｔ）
も、ａｘ＋ｂｙ＝ｇの解になります。
　　　９２４×(－７＋３０ｔ)＋３６０×(１８－７７ｔ)
　　＝－６４６８＋２７７２０ｔ＋６４８０－２７７２０ｔ
　　＝１２

プログラム

　　function ex_euclid(a, b, gcd, x, y){
　　　　r0 = a;　r1 = b;
　　　　x0 = 1;　x1 = 0;
　　　　y0 = 0;　y1 = 1;
　　　　while (r1>0) {
　　　　　　q = Math.floor(r0/r1);
　　　　　　r2 = r0 % r1;
　　　　　　x2 = x0 - q*x1;
　　　　　　y2 = y0 - q*y1;
　　　　　　r0 = r1;　r1 = r2;
　　　　　　x0 = x1;　x1 = x2;
　　　　　　y0 = y1;　y1 = y2;
　　　　}
　　　　gcd = r0;
　　　　x = x0;
　　　　y = y0;
　　}

トレース

ａ＝９２４，ｂ＝３６０としたときの経過を示します。

ｑ　　　ｒ０　　ｒ１　　ｒ２　　　ｘ０　ｘ１　ｘ２　　　ｙ０　ｙ１　ｙ２
２　　９２４　３６０　２０４　　　　１　　０　　１　　　　０　　１　－２
１　　３６０　２０４　１５６　　　　０　　１　－１　　　　１　－２　　３
１　　２０４　１５６　　４８　　　　１　－１　　２　　　－２　　３　－５
３　　１５６　　４８　　１２　　　－１　　２　－７　　　　３　－５　１８
４　　　４８　　１２　　　０　　　　２　－７　３０　　　－５　１８－７７
４　　　１２　　　０　　　０　　　－７　３０　３０　　　１８－７７－７７
　　　　　↓　　　↓　　　　　　　　↓　　↓　　　　　　　↓　　↓
　　　　　ｇ 　　終了　　　　　　　 ｘ　 b/g　　　　　　　ｙ　 -a/g
基本解
　　　ａ　　　　ｘ　　　ｂ　　　　ｙ　　 GCD
　　９２４×（－７）＋３６０×（１８）＝１２
一般解
　　　ｂ／ｇ＝ｘ１＝　３０　　ｘ＝－７＋３０ｔ
　　－ａ／ｇ＝ｙ１＝－７７　　ｙ＝１８－７７ｔ

計算プログラム