お題は、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
ソースはここに置いた。
おもしろかったー。
-
PairingHeap は Chris Okasaki の『Purely Functional Data Structures』で紹介されていた、不思議とパフォーマンスが良いという面白いデータ構造。 ↩