Golangで「エラトステネスの篩」で「2.1秒で百万個」の素数を計算できる「無限シーケンス」を作ってみた

※ タイトルは2.1秒となっていますが、0.96秒を達成しました。詳しくは「追記の章」を参照。

## はじめに
[C#で「エラトステネスの篩」で「2.6秒で百万個」の素数を計算できる「無限シーケンス」を作ってみた - tail-island](https://tail-island.github.io/programming/2016/04/15/lazy-primes.html)を参考に、Golangの<b>チャネルとgoroutine</b>と<b>クロージャ</b>をそれぞれ用いた無限素数ジェネレータを実装してみました。タイトルは元記事のオマージュです。

みなさんご存知のエラトステネスの篩は、自然数のリストから素数の倍数を消していくものだと思います。無限に素数を生成するためには遅延評価などを使って無限に素数の倍数を消していく必要があります。しかし、遅延評価がデフォルトではないGolangでは、自然数リストを有限にしか持つことができないため、決め打ちの値までしか素数を生成できませんでした。そこで、今回実装した無限素数ジェネレータは、素数の奇倍数と素数の2倍の数をペアで管理することで、無限の自然数リストを持たなくてもよい設計としました（型や計算機による制限はありますが）。

ちなみに、元記事では[The Genuine Sieve of Eratosthenes](https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf)という論文を参照しているそうです。

## チャネルを用いた実装
``func NewPrimeGenerator() *PrimeGenerator``でgoroutineを起動し無限素数ジェネレータを作成した後、``func (p *PrimeGenerator) Next() uint32``経由でチャネルから素数を取得しています。

```go:prime_generator_with_channel.go
package main
 
import (
	"fmt"
)
 
type PrimeGenerator struct {
	primeChan chan uint32
}
 
func NewPrimeGenerator() *PrimeGenerator {
	primeGenerator := PrimeGenerator{
		primeChan: make(chan uint32),
	}
 
	go primeGenerator.start()
 
	return &primeGenerator
}
 
func (p *PrimeGenerator) start() {
	// Key: 元の素数の奇倍数, Value: 元の素数の2倍の数
	multiples := map[uint32]uint32{}
	
	// 2は偶素数なので例外扱い
	p.primeChan <- 2

	// 3以上の奇数の中で素数であるものを探す
	for num := uint32(3); ; num += 2 {
		factor, hasFactor := multiples[num]

		// hasFactorがtrueならば、numはfactor/2の倍数なので素数ではない
		if hasFactor {
			delete(multiples, num)
		} else {
			// 後に追加する倍数を奇数に限定するための調節
			factor = num * 2
		}
 
		// 新たな倍数を追加
		for newNum := num + factor; ; newNum += factor {
			_, hasNewFactor := multiples[newNum]

			// hasNewFactorがtrueならば、newNumは既に登録された倍数なのでスルー
			if !hasNewFactor {
				multiples[newNum] = factor
				break
			}
		}
 
		if !hasFactor {
			p.primeChan <- num
		}
	}
}
 
func (p *PrimeGenerator) Next() uint32 {
	return <- p.primeChan
}
 
func main() {
	n := 1000000
	primeGenerator := NewPrimeGenerator()
	for i := 1; i < n; i++ {
		primeGenerator.Next()
	}
	fmt.Printf("%d番目の素数: %d", n, primeGenerator.Next())
}
```

## クロージャを用いた実装
上記のチャネルを用いた実装は、以下の通りクロージャを用いても実装できます。チャネルでの実装をクロージャでの実装に変換できる場合、一般的に後者の方が実行速度は早いそうです。

```go:prime_generator_with_clojure.go
package main
import (
	"fmt"
)
 
type PrimeGenerator struct {
	Next func() uint32
}
 
func NewPrimeGenerator() *PrimeGenerator {
	primeGenerator := PrimeGenerator{}
 
	multiples := map[uint32]uint32{}
	var num uint32
	primeGenerator.Next = func() uint32 {
		if num == 0 {
			num = 1
			return 2
		}
 
		for ;; {
		    num += 2
 
			factor, hasFactor := multiples[num]
			if hasFactor {
				delete(multiples, num)
			} else {
				factor = num * 2
			}
 
			for newNum := num + factor; ; newNum += factor {
				_, hasNewFactor := multiples[newNum]
				if !hasNewFactor {
					multiples[newNum] = factor
					break
				}
			}
 
			if !hasFactor {
				return num
			}
		}
	}
 
	return &primeGenerator
}
 
func main() {
	n := 1000000
	primeGenerator := NewPrimeGenerator()
	for i := 1; i < n; i++ {
		primeGenerator.Next()
	}
	fmt.Printf("%d番目の素数: %d", n, primeGenerator.Next())
}
```

## ベンチマーク

### 実行環境
他の言語との速度比較を行いやすいように、オンラインプログラム実行サービスのIdeoneを使用しました。Ideoneはサポートされている言語の種類が多く、言語のバージョンも比較的新しいので便利ですね。

- Ideone (2017/02/12時点)
- Go (go 1.7.4)

### 結果
実装した無限素数ジェネレータで素数を100万個生成したときの実行時間は、以下の通りになりました。

|| with channel | with clojure |
|:-:| :-: | :-: |
|Time|2.28s|2.1s|
|Site| https://ideone.com/DEfeV2 | https://ideone.com/zjVCgk |

やはり、クロージャーのほうが早いですね。

## まとめ
今回は、チャネルとクロージャを用いて素数無限シーケンスを作成しました。クロージャのほうがチャネルより早い結果となりましたが、チャネルのほうはバッファを使うことでもう少し早くできそうですね。

## 追記の章

[追記1] 2.1秒 → 0.96秒
@antimon2さんの[Golangで「エラトステネスの篩」で「1.3秒で百万個」の素数を計算できる「無限シーケンス」を作ってみた](http://qiita.com/antimon2/items/cada59fb3b1f028f4fb3)を参考に、さらに5, 7の倍数をスキップすることで、100万個の素数を<b>[0.96 sec](https://ideone.com/qBXbRM)</b>で生成することができました。なので、この記事のタイトルの実行時間は高速化前ものになります。

[追記2] 様々な言語で実装
[Extensible prime generator - rosettacode.org](https://rosettacode.org/wiki/Extensible_prime_generator)に[The Genuine Sieve of Eratosthenes](https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf)の実装がいくつか見つかったのでご紹介。

- 実装言語と100万個の素数を生成する時間

|言語|Time|
|:--:|:--|
|[D](https://rosettacode.org/wiki/Extensible_prime_generator#D)|[1.9 sec](https://ideone.com/n5EfqI)|
|[Go](https://rosettacode.org/wiki/Extensible_prime_generator#Go)|[4.54 sec](https://ideone.com/pc4QWQ)|
|[Clojure](https://rosettacode.org/wiki/Extensible_prime_generator#Clojure)|[TLE(5 sec+)](https://ideone.com/NMkH39)|