Help us understand the problem. What is going on with this article?

シクシク素数アドベントカレンダー Scala編

お題は、4か 9を含む素数を指定個数だけ列挙するというものだった。これを下記の趣向で解いてみたい。

  • 素数生成のロジックをあえて再発明してみる
  • ミュータブルなロジックとイミュータブルなロジックの二通りで書く
  • なるべくガチでやってみる

共通コード

共通コードは下記のように書いた。このprimeメソッドを実装するクラスを、二通りに書いてみる。

trait Prime4949Solver {
  def primes: Stream[Int]
  def is4949Prime(n: Int): Boolean = {
    def loop(m: Int): Boolean = (m / 10, m % 10) match {
      case (_, 4) | (_, 9) => true
      case (0, _)          => false
      case (d, _)          => loop(d)
    }
    loop(n)
  }
  def main(args: Array[String]): Unit = {
    val n = args.head.toInt
    println(primes.filter(is4949Prime).take(n).mkString(","))
  }
}

is4949Primeは、文字列操作を使えば短く書けるが、ここでは末尾再帰で実装した。

ミュータブル版: エラトステネスの篩

ミュータブル版は、メモリ領域を書き換えながら合成数を消去していく方式でエラトステネスの篩を実装してみる。

まず篩にかける範囲をざっくり決めたい。100個のシクシク素数が収まるような充分に大きな範囲を適当にとってもおそらく問題ないが、一応見積もってみる。

ある素数がシクシク素数である確率を考えてみると、1桁の素数の場合、確率0なのは自明。n>1桁の場合、一桁目は {1, 3, 7, 9} から 4分の1の確率で9が現れ、二桁目より上の桁にはそれぞれ10分の1の確率で4か9が現れると仮定すると、シクシク素数である確率$p_1$と、そうでない確率$p_2$ は以下の式で求められる。

\begin{pmatrix}
  p_1\\
  p_2
 \end{pmatrix}
=
\begin{pmatrix}
  1 & \frac{2}{10} \\
  0 & \frac{8}{10}
 \end{pmatrix}^{n-1}
\begin{pmatrix}
  \frac{1}{4} \\
  \frac{3}{4}
 \end{pmatrix}

また $10^n$ より小さい素数の個数 $\pi(10^n)$を、素数定理より $10^n / \log 10^n$ で近似すると、下表のように個数が見積れる。

桁数 p1 π 素数 シクシク素数
1 0.0 4.3 約4個 0個
2 0.4 21.7 約17個 約7個
3 0.52 144.8 約123個 約64個
4 0.62 1085.7 約941個 約584個

桁数が小さいのでだいぶ粗い近似になるが、4桁もあれば100個までのシクシク素数を含むのに十分だとわかる。

これを踏まえた範囲で、mutable.BitSet を用いた篩を実装すると以下のようになる。

import scala.collection.mutable

object Main4949Prime extends Prime4949Solver {
  def sieve(limit: Int): Stream[Int] = {
    val bs = mutable.BitSet(limit)
    (2 to limit) foreach (bs += _)
    (2 to math.sqrt(limit).toInt) foreach { n =>
      (n * n to limit by n) foreach (bs -= _)
    }
    bs.toStream
  }
  override def primes = sieve(9999)
}

イミュータブル版: 遅延評価 + Pairing Heap + Wheel Sieve

Haskell の Data.Numbers.Primes では、Pairing Heap1 と Wheel Sieve を併用した素数の無限リストが提供されている。ここではこれを簡略化したものを Scala の Stream を使って実装する。

まず Pairing Heap は以下のように書いた。

sealed trait PairingHeap {
  def merge(another: PairingHeap): PairingHeap
}
case object Empty extends PairingHeap {
  def merge(another: PairingHeap): PairingHeap = another
}
case class Tree(ns: Stream[Int], phs: List[PairingHeap]) extends PairingHeap {
  private def join(t: Tree): PairingHeap =
    if (ns.head <= t.ns.head) Tree(ns, t :: phs) else t.join(this)

  def merge(another: PairingHeap): PairingHeap = another match {
    case Empty   => this
    case t: Tree => join(t)
  }
}
object PairingHeap {
  def enqueue(ns: Stream[Int], q: PairingHeap): PairingHeap = Tree(ns, Nil).merge(q)
  def mergeAll(phs: List[PairingHeap]): PairingHeap = phs match {
    case Nil               => Empty
    case ph1 :: Nil        => ph1
    case ph1 :: ph2 :: pht => ph1 merge ph2 merge mergeAll(pht)
  }
}

Wheel Sieve(Factorization) とは、素数列の初めの既知の小さな部分列から素数候補の数列を作って、効率的に素数を選択する技術。

例えば、{2, 3, 5} なら、2 と 3 から周期 2 * 3 = 6 となるような差分の列 {2, 4}を作り、これを 5以降に繰り返し加算して {2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, ...}といった素数候補列を生成する。(実用的には、最初の {2, 3, 5, 7, 11 ,13, 17} くらいが丁度良いらしいが、ここでは簡単のため {2, 3, 5} にした。)

この素数候補の数列は、下記の要領で合成数の数列を作るのにも使う(各行が遅延リストで実装される)。

{4, 6, 10, 14, 22, 26, ...}
  { 9, 15, 21, 33, 39, ...}
     { 25, 35, 55, 65, ...}
         { 49, 77, 91, ...}
         ...

この「無限数列の無限列」を Pairing Heap に入れて先頭から順に取り出すと、小さい順に並んだ合成数の無限数列が得られる。これを素数候補の数列と突き合わせて取捨選択することで、無限素数リストを生成できる。

以下のようなコードになる。

object WheelSieveSolver extends Prime4949Solver {
  import PairingHeap._

  private def spin(start: Int, wheel: Stream[Int]): Stream[Int] = {
    def repeat: Stream[Int] = wheel append repeat
    repeat.scan(start)(_ + _)
  }
  private def composites(ns: Stream[Int]): Stream[Int] = {
    def loop(ms: Stream[Int]): Stream[Int] = (ns.head * ms.head) #:: loop(ms.tail)
    loop(ns)
  }
  private def sieve(ns: Stream[Int], ph: PairingHeap): Stream[Int] = (ns, ph) match {
    case (n #:: nt, Empty)              => n #:: sieve(nt, enqueue(composites(ns), Empty))
    case (n #:: nt, Tree(m #:: ms, qs)) =>
      if      (n < m) n #:: sieve(nt, enqueue(composites(ns), ph          ))
      else if (n > m)       sieve(ns, enqueue(ms,             mergeAll(qs)))
      else                  sieve(nt, enqueue(ms,             mergeAll(qs)))
  }
  override def primes: Stream[Int] = sieve(Stream(2, 3) append spin(5, Stream(2, 4)), Empty)
}

実行結果

❯❯❯ scala WheelSieveSolver 100
19,29,41,43,47,59,79,89,97,109,139,149,179,191,193,197,199,229,239,241,269,293,347,349,359,379,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,509,541,547,569,593,599,619,641,643,647,659,691,709,719,739,743,769,797,809,829,839,859,907,911,919,929,937,941,947,953,967,971,977,983,991,997,1009,1019,1039,1049,1069,1091,1093,1097,1109,1129,1193,1229,1249,1259,1279,1289,1291,1297,1319
❯❯❯ scala EratosthenesSieveSolver 100
19,29,41,43,47,59,79,89,97,109,139,149,179,191,193,197,199,229,239,241,269,293,347,349,359,379,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,509,541,547,569,593,599,619,641,643,647,659,691,709,719,739,743,769,797,809,829,839,859,907,911,919,929,937,941,947,953,967,971,977,983,991,997,1009,1019,1039,1049,1069,1091,1093,1097,1109,1129,1193,1229,1249,1259,1279,1289,1291,1297,1319

ソースはここに置いた。
おもしろかったー。


  1. PairingHeap は Chris Okasaki の『Purely Functional Data Structures』で紹介されていた、不思議とパフォーマンスが良いという面白いデータ構造。 

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away