素因数分解 - 倭算数理研究所

Scala で『Java によるアルゴリズム事典』のコードを書き換えてみようシリーズ（目次）。　前回、素数列を生成するアルゴリズムを見ましたが、今回は素因数分解 (factorization into primes) をするアルゴリズムを見ていきます。

参考

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

アルゴリズム

ある自然数を素因数分解するには素数で割っていけばいいのですが、素数を生成するのにもコストがかかるので、2と3以上の奇数で順々に割り切れるかどうか（素因数であるかどうか）を試していく方が簡単で効率的です。　例えば、全ての素因数3を抜き出した後では、9で割っても必ず余りが出るので素因数には影響が出ません。　詳しくは実装コード参照。

素数の代わりに２と3以上の奇数で割っていくと数が大きくなるにつれて無駄撃ちが多くなるので、2, 3 以降では 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, ... のように2, 4 ずつ交互に離れた奇数で割っていく方が能率的になります。

また、素数の場合と同じく、素因数の候補は元の自然数の（正の）平方根以下でしかあり得ないので、それを超える数で割り切れるかどうかを試す必要はありません。

実装

では実装コード。　まずは素因数がフラットに並んだ Seq で返される flatFactorize() メソッドを定義しましょう。　「素因数がフラットに並んだ」とは、例えば24を素因数分解したときに Seq(2, 2, 2, 3) のように重複した素因数が重複の個数だけ並んで出力されることを表しています。

  def flatFactorize[I: Integral](n: I): Seq[I] = {

    def factorize(n: I, ps: Stream[I]): Stream[I] =
      n match {
        case 1 => Stream()
        case _ => ps.head match {  // ps.head が次の素因数候補
          case p if p * p > n  => Stream(n)  // p 以上 n 未満の素因数なし => n 自体が素因数
          case p if n % p == 0 => p #:: factorize(n /~ p, ps)  // p が素因数
          case _               => factorize(n, ps.tail)  // p は素因数ではない => 次の候補へ
        }
      }

    // 素因数候補の Stream
    val ps = 2 #:: IntegralSequence.from[I](3, 2)

    factorize(n, ps)
  }

n が1なら空の Stream を返します（1は素数ではないので）。
#:: は Stream トレイトのメソッド（「:」で終わるので右結合）で、引数 (p) を先頭にした遅延評価の Stream を生成します。
/~ は spire での整数の割り算です。　通常の Scala や Java での整数の / と同じです。
ps は素因数の候補の列です。　ここでは2の後に3以上の奇数が続く Stream です。
IntegralSequence.from() メソッドは素数の記事で定義したユーティリティ・メソッドで、今の場合、3で始まり2ずつ増える等差数列を生成します。　標準 API の Stream.from() メソッドと同じですが、spire の Integral に対しても動くというだけです。

もう少し効率的にするには、素数候補の ps を以下の定義で変えるだけ OK：

    val ps = 2 #:: 3 #:: IntegralSequence.from[I](5, 6).flatMap(i => Seq(i, i+2))

以上で素因数分解を行うコードは終わりですが、結果の Seq がちょっと扱いにくいので、少々イジって (素因数, 重複数) というタプルの Stream に変換しておきましょう。　素因数は小さいもの順に並んでいるのでそんなに難しくなく変換できます：

  def groupSequentialDuplicates[A](seq: Seq[A]): Stream[(A, Int)] =
    seq match {
      case Nil => Stream()
      case head +: _ =>
        val dup = seq.takeWhile(_ == head)
        val n = dup.length
        (head, n) #:: groupSequentialDuplicates(seq.drop(n))
    }

この groupSequentialDuplicates() メソッド（長い名前だけど直接は使わない前提で名付けてます）と先ほどの flatFactorize() メソッドを使って、factorize() メソッドを以下のように定義しましょう：

  def factorize[I: Integral](n: I): Seq[(I, Int)] =
    groupSequentialDuplicates(flatFactorize(n))

これで素因数分解のコード完成。

サンプル・コード

では実際に使って見ましょう。　一応 flatFactorize() も試しておきます：

// flatFactorize() メソッド
assert( flatFactorize(1) = Seq() )  // 1は素因数なし
assert( flatFactorize(2*3*4*5) == Seq(2, 2, 2, 3, 5) )
assert( flatFactorize(10*11*12*13*14*15) == Seq(2, 2, 2, 2, 3, 3, 5, 5, 7, 11, 13) ) )
assert( flatFactorize(1024) == Seq.fill(10)(2) )

// factorize() メソッド
assert( factorize(1) = Seq() )
assert( factorize(2*3*4*5) == Seq((2, 3), (3, 1), (5, 1)) )  // 2*3*4*5 = 2^3 * 3^1 * 5^1
assert( factorize(10*11*12*13*14*15) == Seq((2, 4), (3, 2), (5, 2), (7, 1), (11, 1), (13, 1)) )
assert( factorize(1024) == Seq((2, 10)) )

大丈夫そうですかね。　ジェネリックなコードになってるのでもちろん BigInt でも動きます：

import scala.language.postfixOps
import spire.implicits._  // BigInt に対して「!」で階乗を計算できるようにする

// 素因数分解の結果を見やすく整形
def s[A](seq: Seq[(A, Int)]): String = seq.map{
  case (p, 1) => p.toString
  case (p, e) => p + "^" + e
}.mkString(" * ")

// 10! の素因数分解
assert( s(factorize(BigInt(10)!)) == "2^8 * 3^4 * 5^2 * 7" )

// 100! の素因数分解ぃ！
assert( s(factorize(BigInt(100)!)) == 
  "2^97 * 3^48 * 5^24 * 7^16 * 11^9 * 13^7 * 17^5 * 19^5 * 23^4 * 29^3 * 31^3 * 37^2 * 41^2 * 43^2 * 47^2 * 53 * 59 * 61 * 67 * 71 * 73 * 79 * 83 * 89 * 97" )