読者です 読者をやめる 読者になる 読者になる

倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

素数の Lucas テスト

Scalaアルゴリズム事典 数論

今回はメルセンヌ数(Mersenne Number)が素数であるかどうかを判定するリュカ・テスト (Lucas–Lehmer primality test) を見ていきます(目次)。

メルセンヌ数とは  { p }素数として

  { \displaystyle\begin{align*}
  M_p = 2^p - 1
\end{align*}}

の形の数のことです。 ただし、この記事では  { p }素数でなくても  { M_p }メルセンヌ数ということにします。

参考

アルゴリズム

リュカ・テストは、奇素数  { p } に対して

  { \displaystyle\begin{align*}
  \begin{cases}
    x_0 \;\!\! &= 4 \\
    x_{i+1} \!\!&= (x_i^2 -2) \mod M_p \qquad (i \geqq 0)
  \end{cases}
\end{align*}}

で定義される数列  { x_i } に対して

  { \displaystyle\begin{align*}
  x_{p-2} = 0    &\quad\Longrightarrow\quad M_p\,\textrm{は素数} \\
  x_{p-2} \ne 0 &\quad\Longrightarrow\quad M_p\,\textrm{は合成数}
\end{align*}}

が成り立つことを使って、メルセンヌ数  { M_p }素数かどうかを判定します。 証明は 『wikipedia:リュカ–レーマー・テストの証明』 に日本語で書いてくれているので、興味のある方はどうぞ。

補足
リュカ・テストが成り立つ証明は上記のリンク先に載ってますが、その記事では  { x_i } の代わりに以下で定義される  { S_i } を使っています:

  { \displaystyle\begin{align*}
  \begin{cases}
    S_0 \;\!\! &= 4 \\
    S_{i+1} \!\!&= S_i^2 -2 \qquad (i \geqq 0)
  \end{cases}
\end{align*}}

で定義される数列  { S_i } に対して

  { \displaystyle\begin{align*}
  S_{p-2} \equiv 0 \mod M_p &\quad\Longrightarrow\quad M_p\,\textrm{は素数} \\
  S_{p-2} \not\equiv 0 \mod M_p &\quad\Longrightarrow\quad M_p\,\textrm{は合成数}
\end{align*}}

 { \mod M_p } を漸化式の各評価ごとにとるか、最後にとるかの違いですが、結果には影響しません。 プログラム中では漸化式の途中でとっていかないと数列の各項が大きくなりすぎるので、 { x_i } の方の定義を用います。

また、リンク先の記事の「必要性」の箇所に

 { p } が奇素数であり、かつ  { Q }素数 { S_{p-2} \equiv 0 \mod Q }

( { Q = M_p }) とありますが、対偶をとると

 { S_{p-2} \not\equiv 0 \mod Q } { p } が奇素数でない、または  { Q }素数でない

となります。 これを見ると、もし  { p } が2(偶素数)もしくは合成数なら  { Q }素数でなくてもよさそうに見えます。 しかし、 { p }合成数なら  { p = mn\quad (m,\,n \geqq 2) } とおいて

  { \displaystyle\begin{align*}
  M_p
    &= 2^{mn} - 1 \\
    &= (2^n-1)\left(2^{(m-1)n} + 2^{(m-2)n} + 2^{(m-3)n} + \cdots + 2^{n} + 1\right)
\end{align*}}

と変形でき、どちらの因子も3以上なので素因数分解できることが分かります。 つまり  { p }合成数なら  { Q }合成数となります。  { p } が2の場合は  { S_{p-2} = S_0 = 4 \ne 0 } よりテストが通りませんが2は素数なので、この場合はリュカ・テストとは別に素数判定する必要があります。

 { S_i } の一般項
リュカ・テストの証明はちょっと込み入ってるようですが、 { S_i } の一般項はそんなに難しくなく出せるのでちょっとやってみましょう。 数列  { \omega_i }

  { \displaystyle\begin{align*}
  S_i = \omega_i + \frac{1}{\omega_i}
\end{align*}}

で定義すると、 { S_i } の漸化式から

  { \displaystyle\begin{align*}
  S_{i+1}
    &= S_i^2 - 2 \\
    &= \left(\omega_i + \frac{1}{\omega_i}\right)^2 - 2 \\
    &= \omega_i^2 + \frac{1}{\omega_i^2}
\end{align*}}

となります。 これと  { S_{i+1} = \omega_{i+1} + \dfrac{1}{\omega_{i+1}} } を見比べると、 { \omega_i } の漸化式を

  { \displaystyle\begin{align*}
  \omega_{i+1} = \omega_i^2
\end{align*}}

と定義すればよいことが分かります。 また、初項  { \omega_0 }

  { \displaystyle\begin{align*}
  &4 = S_0 = \omega_0 + \frac{1}{\omega_0} \\
  &\omega_0^2 - 4\omega_0 + 1 = 0 \\
  \therefore \, &\omega_0 = 2 \pm \sqrt{3}
\end{align*}}

複号のうちどちらをとっても  { S_i } には影響しないので  { \omega_0 = 2 + \sqrt{3} } としましょう。 このとき

  { \displaystyle\begin{align*}
  \omega_i &= (2+\sqrt{3})^{2^i} \\[2mm]
  \therefore\, S_i &= (2+\sqrt{3})^{2^i} + (2-\sqrt{3})^{2^i}
\end{align*}}

となって、 { S_i } の一般項が求まりました*1

実装1

では実装してみましょう。 メルセンヌ数はそんなに大きくない整数  { p } に対してもすぐに Int や Long の最大値を超えるので(それぞれ 31, 63 まで)BigInt を使うようにしましょう:

import spire.implicits._

case class MersenneNumber(p: Int){
  require(p >= 0)

  lazy val toBigInt: BigInt = BigInt(2)**p - 1

  def isPrime: Boolean = p match {
    case 0 | 1 => false
    case 2 => true  // p = 2 の場合
    case _ =>
      val m = toBigInt
      val seq = Stream.iterate[BigInt](4)(x => (x * x - 2) % m)
      seq(p - 2) == 0
  }
}
  • 初項と漸化式が与えられている数列の Seq を生成したい場合、Stream.iterate メソッドを使うのがいいでしょう。
  • 数学の mod は、コード中では言わずと知れた「%」で。

まぁ別にクラスにする必要はなかったかも知れませんが。 使い方は以下のようになります:

// p = 3 の場合
val mn3 = MersenneNumber(3)
assert( mn3.toBigInt == BigInt(7) )
assert( mn3.isPrime == true )

// p = 11 の場合(メルセンヌ数は素数ではない)
val mn11 = MersenneNumber(11)
assert( mn11.toBigInt == BigInt(2047) )
assert( mn11.isPrime == false )

// p が 640 まででメルセンヌ数が素数になるもの
assert( (0 to 640).filter(MersenneNumber(_).isPrime) ==
          Seq(2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607) )

p = 607 のときのメルセンヌ数ってかなり大きい数になるでしょうけど、このコードは大した時間かからずに処理が帰ってきます。 おぉすごい!

実装2

Java によるアルゴリズム事典』には、BigInt (BigInteger) を使った実装以外に、配列を使った実装も載ってあるのですが、そのアルゴリズムScala の Seq 等で書いてみるとかなりパフォーマンスがクソなことになったので、ここではメルセンヌ数を法とする和・差・積を定義した MersenneModuloUInt というクラスを作って、計算をどのようにするかの簡単な例としたいと思います。

ぶっちゃけ、BigInt の実装で何の問題もないので、これ以降は特に読む必要はないかと思います。

メルセンヌ数を法とする代数
メルセンヌ数  { M_p = 2^p-1 } { p }素数でなくてもよい)を法とした代数の要素は0から  { 2^p-2 } までの  { 2^p-1 } 個の数で表せるので、 { p } ビットで表すことができます。  { p } ビットでは  { 2^p } 個の数を表せるので、実際には1つ余計で(右下の添字2で2進数を表し、合同  { \equiv } も = を使うことにして)

  { \displaystyle\begin{align*}
  \underbrace{111\cdots1}_{p\textrm{個}}{}_2 = \underbrace{000\cdots0}_{p\textrm{個}}{}_2 = 0
\end{align*}}

を同一視します。

最大の数  { 2^p-2 }

  { \displaystyle\begin{align*}
  2^p-2 = \underbrace{111\cdots1}_{(p-1)\textrm{個}}0_2
\end{align*}}

となります。

 { M_p } を法とした場合、ビット数が  { p } より大きい数を mod (プログラムでは %)によって  { p } ビットに収まる数に変換するのが、比較的簡単なビット演算で計算できます。 上で見たように  { p } 個のビットが全て1ならその数は0となるので、その次の数を考えて

  { \displaystyle\begin{align*}
  1\underbrace{000\cdots0}_{p\textrm{個}}{}_2 &= \underbrace{000\cdots0}_{(p-1)\textrm{個}}{1}_2 = 1
\end{align*}}

となりますが、これは下位  { p } 個の0を取り除くことに対応しています。 ビット演算で言えば、符号なしで  { p } だけ右シフトすることに相当します。

負の数(というか和に関する逆元)は単に全ビットを反転すればいいだけです。 実際、全ビットが逆の2数を足せば全ビットが1の数になりますが、これは0に等しいんでしたね。 ただし、0の場合はそのまま0です。 全ビットを反転して全て1になってもそれは0です。

ちょっと具体例
もう少し理解を深めるために、具体的な数字で計算してみましょう。  { p = 3 } のとき  { M_3 = 7 } ですね。 このとき、最大の数は  { 110_2 = 6 }、また  { 111_2 = 7 \equiv 0 \!\!\mod 7 } となります。  { p } ビットより大きい数に対しては

  { \displaystyle\begin{align*}
  100110_2
    &= 100\underline{000}_2 + 110_2 \\
    &= 100_2 + 110_2 \\
    &= 1010_2 \\
    &= 1\underline{000}_2 + 10_2 \\
    &= 1_2 + 10_2 \\
    &= 11_2
\end{align*}}

 { 100110_2 = 38 \equiv 3 \!\!\mod 7 } なので確かに合ってますね。 上記の計算では少々助長に途中式を書いてますが、要は右から  { p } 個ごとにビットを区切って足せば OK です:

  { \displaystyle\begin{align*}
  100110_2
    &= 100_2 + 110_2 \\
    &= 1010_2 \\
    &= 1_2 + 10_2 \\
    &= 11_2
\end{align*}}

負の数は、例えば  { 101_2 (= 5) } に対して、全ビットを反転した  { 010_2 = 2 } となります。 実際にこの2数を加えてみると

  { \displaystyle\begin{align*}
  101_2 + 010_2
    &= 111_2 \\
    &= 000_2
\end{align*}}

となって、確かに和に関して逆元になっています。

コード
せっかく代数の軽い説明をしたので、その代数の演算が定義されたクラスを作りましょう。 何桁でも扱えるようにするのは非常に大変なので、Int で表せる32ビットまでを扱うことにします( { p \leqq 32 })。 実質的に非負の整数しか出てこないので、プロパティとしては Int ではなく spire に定義されている UInt を使う方がよさそうです( { p = 32 } も扱えるので)。 そのクラスを MersenneModularUInt として各種演算を定義しますが、まずはこのオブジェクトを生成する MersenneModulus クラスを定義しましょう:

import spire.math.UInt
import spire.implicits._
import scala.annotation.tailrec

case class MersenneModulus(p: Int){
  require(2 <= p && p <= 32)

  val mask64 = 2L**p-1L
  val mask32 = UInt(mask64.toInt)
    // For p = 3 :
    //  mask = 00...0111
    // ~mask = 11...1000

  /** UInt から MersenneModularUInt オブジェクトを生成するメソッド */
  def fromUInt(n: UInt): MersenneModularUInt = {
    @tailrec
    def mod(n: UInt): UInt = n & ~mask32 match {
      case x if x == UInt(0) =>
        if(n == mask32) UInt(0) else n
      case x => mod((n & mask32) + (n >>> p))
    }

    MersenneModularUInt(mod(n), this)
  }

  /** Int から MersenneModularUInt オブジェクトを生成する apply ファクトリ・メソッド */
  def apply(n: Int): MersenneModularUInt =
    if(n < 0) -fromUInt(UInt(-n))
    else       fromUInt(UInt(n))

  /** Long から MersenneModularUInt オブジェクトを生成する apply ファクトリ・メソッド */
  def apply(n: Long): MersenneModularUInt = n match {
    case _ if n < 0L => -apply(-n)
    case _ =>
      @tailrec
      def mod(n: Long): Long = n & ~mask64 match {
        case 0L =>
          if(n == mask64) 0L else n
        case x => mod((n & mask64) + (n >>> p))
      }

      MersenneModularUInt(UInt(mod(n)), this)
  }
}
  • Int 値や Long 値をとって、後で定義する MersenneModularUInt オブジェクトを生成する apply() ファクトリ・メソッドを定義しています。
  • UInt 値をとる apply() ファクトリ・メソッドも定義しようとしたんですが、Int をとるものが定義されているとダメなようなので、fromUInt() メソッドとして定義しています。
  • p ビットを超える数や負の数を、ビットシフトを使って p ビットで表せる数に変換しています。

さて、次は各種演算が定義された MersenneModularUInt クラス:

import spire.math.UInt
import spire.implicits._

case class MersenneModularUInt(digit: UInt, modulus: MersenneModulus){

  def p: Int = modulus.p

  def isZero: Boolean = this.digit == UInt(0)
  def isOne : Boolean = this.digit == UInt(1)

  def unary_- : MersenneModularUInt =
    if(this.isZero) this
    else            MersenneModularUInt(modulus.mask32 & ~digit, modulus)

  def +(that: MersenneModularUInt): MersenneModularUInt =
    if(this.isZero) that
    else if(that.isZero) this
    else calculateSum(that)

  /** 足し算実行 */
  private def calculateSum(that: MersenneModularUInt): MersenneModularUInt =
    this.modulus(this.digit.toLong + that.digit.toLong)

  def -(that: MersenneModularUInt): MersenneModularUInt =
    if(this.isZero) -that
    else if(that.isZero) this
    else this + (-that)

  def *(that: MersenneModularUInt): MersenneModularUInt =
    if(this.isZero | that.isOne) this
    else if(this.isOne | that.isZero) that
    else calculateProduct(that)

  /** 掛け算実行 */
  private def calculateProduct(that: MersenneModularUInt): MersenneModularUInt =
    this.modulus(this.digit.toLong * that.digit.toLong)

  def toInt: Int = this.digit.toInt
  def toLong: Long = this.digit.toLong
  def toBigInt: BigInt = this.digit.toBigInt
  override def toString: String = s"$digit (mod 2^$p-1)"
}
  • 演算に0や1が現れた場合には計算せずに済むことが多いので、いろいろ条件分岐していますが、和や積の実質的な計算は calculateXxxx という名前の private メソッドで行っています
  • 数値として保持しているのは UInt ですが、和や積の計算途中では結果の桁が溢れると困るので Long に変換してから計算しています。
  • p ビット以下への変換は先に定義した MersenneModulus クラスの apply() ファクトリ・メソッドに任せてます。

さて、この MersenneModularUInt クラスを使ってリュカ・テストを実装すると以下のようになります(まぁ、ほとんど BigInt を使った実装と同じですが):

object UIntMersenneNumber{

  def isPrime(p: Int): Boolean = p match {
    case 0 | 1 => false
    case 2 => true
    case _ =>
      val mmod = MersenneModulus(p)
      val v2 = mmod(2)
      val xs = Stream.iterate(mmod(4))(x => x*x - v2)
      xs(p-2).isZero
  }
}

この isPrime メソッドを使って、32までのメルセンヌ素数を全て求めてみると

assert( (0 to 32).filter(UIntMersenneNumber.isPrime(_)) == Seq(2, 3, 5, 7, 13, 17, 19, 31) )

となってキチンと答えが出ました。

ここで実装した MersenneModularUInt クラスでは 32 以下の p しか扱えませんが、spire の Natural の実装などを参考にすると任意の p に対して計算できるようにすることもできるかと思います。 ただ、Natural のソースを見ると「パフォーマンスを気にする場合は BigInt を使え」とかコメントされてるので、あんまり作る気が起きてきませんがw

【修正】

  • BigInt の実装でのサンプルコード部分が間違っていたので修正しました。
  • 0の和に関する逆元はそのままであることを追記しました。
  • MersenneModulo を MersenneModulus に、MersenneModuloUInt を MersenneModularUInt に名前を変更しました。

Scalaスケーラブルプログラミング第3版

Scalaスケーラブルプログラミング第3版

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

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

*1:どちらかと言うと、この一般項がありきで、コードに落とし込むために漸化式に書き換えたのではないかという気がしますがどうなんでしょうかね。 というか、証明では漸化式ではなく一般項を使ってるのでそういうことなんでしょう。

素因数分解

Scalaアルゴリズム事典 数論

前回素数列を生成するアルゴリズムを見ましたが、今回は素因数分解 (factorization into primes) をするアルゴリズムを見ていきます(目次)。

参考

アルゴリズム

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

素数の代わりに2と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 での整数の割り算です。 通常の ScalaJava での整数の / と同じです。
  • 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, 3, 5) )
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" )

合ってそうかな?

【修正】

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

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

素数

Scalaアルゴリズム事典 数論 数列

ちょっと別用素数 (prime number) を生成するコードを書いたので、こちらのブログで整理しておきます(目次)。 別用記事で見たように Int や Long に対してなら2行で書けますが、ここでは spire の Integral を使ってジェネリックなコードにしましょう。

アルゴリズム

ある自然数が与えられたとき、それが素数かどうかを判定するには、その自然数より小さい全ての素数で割りきれないかどうかを見ます。 実際には、その自然数の(正の)平方根以下の素数で割りきれなければ大丈夫です。 なぜなら、平方根以上の素数で割ったときに割り切れるなら、その商は平方根以下で、その素因数は平方根以下の素数になるからです。

実装

準備
まずはちょっとユーティリティメソッドを定義しておきましょう。 標準 API の Stream.from() メソッドを使えば、無限に続く等差数列を与える Stream オブジェクトを簡単に作れますが、spire の Integral 型に対して同様のメソッドがないので自作しておきます*1

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

object IntegralSequence{
  def from[I: Integral](start: I): Stream[I] = from(start, 1)
  def from[I: Integral](start: I, step: I): Stream[I] = start #:: from(start + step)
}
  • #:: は Stream トレイトに定義されているメソッドで(「:」で終わるので右結合)、遅延評価で無限 Stream を生成します。 Stream を再帰的に定義するときなどに使います。 詳しくは Stream の Scaladoc を参照のこと。

使い方は Stream.from() メソッドと同じです。 以下では引数を2つとる方しか使いませんが、なんとなく引数1つのもの(公差1の等差数列を生成)も定義しています。

素数を生成するコード
では、本題の素数を生成するコード。 以下では4つのメソッドを定義していますが、素数の Stream を生成するアルゴリズムに関係するのは後の2つ。 ただ、簡単に使うには前の2つのメソッドの方が便利かと思います:

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

object PrimeNumber{

  /** n番目の素数を返す(素数2は0番目) */
  def apply[I: Integral](n: I): I = streamOf.apply(n.toInt)

  /** Int 型の素数列を返す */
  def stream: Stream[Int] = streamOf[Int]

  /** I 型の素数列を返す */
  def streamOf[I: Integral]: Stream[I] =
    2 #:: IntegralSequence.from[I](3, 2).filter(isPrime(_))

  /** 素数かどうかを判定する */
  def isPrime[I: Integral](n: I): Boolean =
    n > 1 && streamOf.takeWhile(p => p*p <= n).forall(n % _ != 0)
}
  • apply() メソッドは引数で指定された番号の素数を返します。 素数2は0番目です。
  • stream メソッドは Int の Stream を返します。
  • streamOf メソッドは型パラメータに指定した型の Stream を返します。 たとえば Long 型の Stream を返す場合は streamOf[Long] とします。
    • IntegralSequence.from() メソッドは「準備」のところで定義したメソッドで、ここでは3以上の奇数(3で始まり2ずつ増える)を返す Stream を生成します。 偶数は2以外に素数がないので試すだけ無駄ですね。
    • #:: メソッドは準備のところでも出てきた無限 Stream を生成するメソッドで、ここでは一見再帰になってませんが、isPrime() が streamOf() を呼び出しているので遅延評価にする必要があります。
    • このコードでは、最小で唯一偶数の素数2を先頭にして、その後に奇数から isPrime によって素数をフィルタリングして返しています。
  • isPrime() メソッドは引数の整数が素数かどうかを判定します。 streamOf で返される素数列のうち、引数の平方根より小さいものだけを取り出して、それらすべての素数が引数を割り切らないかどうかを判定しています。

streamOf も isPrime も実質1行ずつで書けてるので、2行で書けているといってもいいでしょう(棒) ちなみに、Int 型に関して2行で書くとこんな感じです:

  val primes: Stream[Int] = 2 #:: Stream.from(3, 2).filter(isPrime)
  def isPrime(n: Int): Boolean = n > 1 && primes.takeWhile(p => p*p <= n).forall(n % _ != 0)

サンプル・コード

では、上記の(ジェネリックな方の)コードを動かしてみましょう。

// apply() メソッドで n 番目の素数を取得する
assert( PrimeNumber(0) == 2 )
assert( PrimeNumber(10) == 31 )

// 素数の Stream[Int] を返す
assert( PrimeNumber.stream.take(10) == Seq(2, 3, 5, 7, 11, 13, 17, 19, 23, 29) )

// 素数の Stream[Long] を返す
assert( PrimeNumber.streamOf[Long].take(10) == 
  Seq(2L, 3L, 5L, 7L, 11L, 13L, 17L, 19L, 23L, 29L) )

streamOf は型パラメータを渡せば BigInt などでも動きますヨ。

【修正】

  • 当初 IntegerSequence としていたコンパニオン・オブジェクトを IntegralSequence に変更しました。

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

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

*1:別に Stream.from() で生成される Stream に map を施しても