最大公約数 - 倭算数理研究所

Scala で『Java によるアルゴリズム事典』のコードを書き換えてみようシリーズ（目次）。　今回は最大公約数（GCD, Greatest Common Divsor）を求めるアルゴリズムを見ていきます。　最大公約数はアルゴリズムと言えばのユークリッド互除法で求められます。

この記事の内容

この記事の内容
参考
2数の最大公約数
- アルゴリズム
- 実装
3個以上の数の最大公約数
- アルゴリズム
- 実装

参考

『Javaによるアルゴリズム事典』

2数の最大公約数

まずは2個の数字の最大公約数を求めるアルゴリズムから見ていきましょう。

アルゴリズム

${ a,\,b }$ を ${ a \geqq b }$ を満たす正の自然数とします。　 ${ a }$ を ${ b }$ で割った商が ${ q }$ 、余りが ${ r }$ であるとき、

　　 ${ \displaystyle\begin{align*} a = bq + r \qquad (q \geqq 1,\quad 0 \leqq r < b) \end{align*}}$

が成り立ちます。　ここで ${ r \ne 0 }$ のとき、 ${ a,\,b }$ の最大公約数を考えると、上式の左辺 ${ a }$ と右辺第1項 ${ bq }$ の共通の約数はまた右辺第2項 ${ r }$ も約数として持たなければならないので、この最大公約数は ${ b,\,r }$ の最大公約数と等しくなります。　また、 ${ r = 0 }$ のとき、 ${ b }$ 自体が ${ a,\,b }$ の最大公約数になります。

これを踏まえて、2つの正の自然数 ${ a,\,b }$ の最大公約数を求めるためには

${ a }$ を ${ b }$ で割り
余りが0なら、 ${ b }$ が最大公約数
余りが0でないなら、 ${ b }$ を新たな ${ a }$ 、余りを新たな ${ b }$ として手順を繰り返す

とすればいいことがわかります。

例えば2数 91, 35 の最大公約数を求めるには以下の様にします：

　　 ${ \displaystyle\begin{align*} 91 &= \underline{35} \times 2 + \underline{21} \\ 35 &= \underline{21} \times 1 + \underline{14} \\ 21 &= \underline{14} \times 1 + \underline{7} \\ 14 &= \textbf{7} \times 2 \end{align*}}$

よって、91, 35 の最大公約数は7であることが分かります。　今の場合は暗算で結果が分かりますが、素因数分解が大変な大きい2数に対しても、この方法なら機械的な計算で最大公約数を求めることができます。

2つの数が負の整数もとりうる場合、2数の最大公約数はそれらの数の絶対値に対して上記の手順で求めたものになります。

2つの数のうち一方が0の場合、0の約数が0以外のすべての整数である*1（a の約数とは、コード的に「a % n == 0 を満たす整数 n」が定義なので、a = 0 なら n は0以外の全ての整数となる）ので、0でない整数 m と0との最大公約数は m となります。　0 同士の最大公約数は数学的には定義されないみたいですが、これを 0 と定義すると以下の実装が簡単になるのでプログラミングでは通常そうするようです。

実装

では実装 ... の前に、spire では spire.math.Integral 型に最大公約数を求める gcd メソッドが定義されています：

import spire.implicits._

assert( 91.gcd(35) == 7  )
assert( (91 gcd 35) == 7 )

assert( (0 gcd 7) == 7 )
assert( (7 gcd 0) == 7 )
assert( (0 gcd 0) == 0 )

もちろん、どちらかもしくは両方の整数が0でもきちんと動きます。

まぁ、これはこれとして、上記のアルゴリズムを実装してみましょう。　名前は gcd() としましょう。　自然数（0を含む）のみを扱う場合は簡単です：

import spire.math.Integral
import spire.implicits._

@tailrec
def gcd[I: Integral](m: I, n: I): I = n match {
  case 0 => m
  case _ => gcd(n, m % n)
}

引数が負の整数もとりうる場合には、引数の絶対値に対して先ほどのアルゴリズムを適用すればいいだけです。　spire では Integral オブジェクトに対して abs メソッドを呼び出せば絶対値を取得できます：

def gcd[I: Integral](m: I, n: I): I = {
  @tailrec
  def gcdPositive(m: I, n:I): I = n match {
    case 0 => m
    case _ => gcdPositive(n, m % n)
  }

  gcdPositive(m.abs, n.abs)
}

使い方は単なるメソッド呼び出し：

assert( gcd( 91,  35) == 7 )
assert( gcd( 91, -35) == 7 )
assert( gcd(-91,  35) == 7 )
assert( gcd(-91, -35) == 7 )

assert( gcd(0, 7) == 7 )
assert( gcd(7, 0) == 7 )
assert( gcd(0, 0) == 0 )

ちなみに、アルゴリズムの説明で ${ a \geqq b }$ という条件を入れてましたが、上記のコードで m < n だと「m % n == m」となって単に1ステップを費やして引数を入れ替えるだけになるので結果に影響しません。　アルゴリズムの説明の際には、単に ${ q = 0 }$ の場合を書くのが面倒にだったのでこういう条件を入れていただけです。

『Java によるアルゴリズム事典』にはビット演算で効率的なアルゴリズムが書かれてたので、ちょいと Scala & spire で書いてみました：

def gcd1[I: Integral](m: I, n: I): I = {
  require(m > 0 && n > 0)
  if(m == 1 || n == 1) return 1

    @tailrec
    def gcd1(m: I, n: I, d: I): I = (m % 2, n % 2) match {
      case (0, 0) => gcd1(m /~ 2, n /~ 2, d * 2)
      case _ =>
        @tailrec
        def gcd1_(u: I, v: I, t: I): I = t match {
          case 0 => u * d
          case _ =>
            @tailrec
            def dropPrimeFactor2(t: I): I = t.abs % 2 match {
              case 0 => dropPrimeFactor2(t /~ 2)
              case _ => t
            }

            dropPrimeFactor2(t) match {
              case s if s > 0 => gcd1_(s,  v, s - v)
              case s          => gcd1_(u, -s, u + s)
            }
        }

        val t = if(m % 2 == 0) m else -n
        gcd1_(m, n, t)
    }

    gcd1(m, n, 1)
  }

「/~」演算子は、spire.math.Integral での整数に対する除算（Java, Scala では単に /）です。　まぁ、興味ある人はコードの内容を解析してくださいまし。

3個以上の数の最大公約数

次は3個以上の整数に関する最大公約数を求めるアルゴリズムを見ていきましょう。　これは上記の2数に関する最大公約数を求めるアルゴリズムを繰り返し使用して求めます。

アルゴリズム

3個以上の数に対しては、まず2個の整数を選んで前節の方法で最大公約数を求め、次にこの数と3個の中の別の数との最大公約数を求め・・・これを繰り返して最後に残った数が求める最大公約数となります。　ただし、途中で1となった場合、それ以上計算する必要がなく最大公約数は1となります。

実装

では実装。　メソッド名は nGCD() とし、Seq[I] を引数にとるものと可変長引数のものを定義しておきましょう：

import spire.math.Integral
import spire.implicits._

def nGCD[I: Integral](seq: Seq[I]): I = {
  require(seq.nonEmpty)
  nGCD(seq.head, seq.tail:_*)
}

def nGCD[I: Integral](m: I, ns: I*): I = m match {
  case 1 => m  // m はそこまでの最大公約数なので、1なら全体の最大公約数も1
  case _ =>
    ns match {
      case Nil => m.abs  // abs は引数が1つしか渡されなかった場合に絶対値を返すため
      case n +: rest => nGCD(gcd(m, n), rest:_*)
    }
}

なんか再帰で書いてくれと言わんばかりのアルゴリズムですねぇ。　使い方は特に問題ないかと：

assert( gcd(Seq(91, 35, 21)) == 7 )  // Seq 引数
assert( gcd(91, 35, 21) == 7 )  // 可変長引数

assert( gcd(0, 7, 14, 21) == 7 )
assert( gcd( 7) == 7 )
assert( gcd(-7) == 7 )
assert( gcd( 0) == 0 )