エラトステネスのふるい

とは？

エラストテネスが考案した速い素数判定方法。準備が面倒なうえに、最大どこまで判定したいかも決めておかなければならなかったりするが、準備を終えると O(1) で判定できる。

ふるった配列の作り方

調べたい値の最大値を 100 とすれば 101 のサイズで用意し、

• 2 の二倍の 4 から 2 ごとに書き込む
• 3 の二倍の 6 から 3 ごとに書き込む

といった感じで 2, 3, 4, ... と繰り返す。

どこまで繰り返す？

√最大値 つまり √100 だから 10 まででよい。これは 10 を超える数の倍数は 10 以下の素数の倍数になっているから。たとえば 11 の二倍の 22 は、2 の倍数になっている。

ベース値があれば飛ばす (重要)

たとえばベースが 4 のとき 4 にはすでに 2 が書き込まれているので、その先も書き込む必要がないのは明らかである。だから飛ばしていい。

にもかかわらず当初、このスキップ処理を忘れてあんまり速くないアルゴリズムだと勘違いしていた。これはふるった印として false を書き込んでいたのがよくなかった。

ベース値を書き込む

初期値を true として false を書き込むのが削除のイメージにあっているのかもしれないが、初期値を空 (nil) としてベース値を書き込んだほうがわかりやすい。

すでに値があれば書き込まない

上書きした方が速そうだが、上書きしない方がどのベース値のときに書き込んだのかわかりやすい。とにかく最初はわかりやすく書く。

使い方

調べたい値をふるった配列の添字にして空なら素数とみなす。

2 以上にしか対応していない

0 と 1 は素数ではないのに、ふるった配列を 0 や 1 の添字で参照してしまうと空なので素数とみなされてしまう。このような場合は 0 と 1 は決め打ちで「素数ではない」としてしまうか、ふるった配列の 0 と 1 を空でないようにしておくのがよさそう。

シミュレーション

最大 100 まで調べたいので 100 + 1 の大きさの配列を用意する。

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100

2 の二倍の 4 から 2 ごとに書き込む。

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100
2					2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2

3 の二倍の 6 から 3 ごとに書き込む。

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100
3					2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2

4 には 2 があるので飛ばす。4 の二倍の 8 から 4 ずつ書き込んでもいいが空はどこにもないので無駄である。

5 の二倍の 10 から 5 ごとに書き込む。

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100
5					2		2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2

6 には 2 があるので飛ばす。

7 の二倍の 14 から 7 ごとに書き込む。

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100
7					2		2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2		2	7	2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2	7	2		2	3	2		2	5	2	3	2		2	7	2	3	2	5	2		2	3	2

ここまで来ると新しく埋められたのは 49 77 91 しかない。

8 には 2 があるので飛ばす。

9 には 3 があるので飛ばす。

10 には 2 があるので飛ばす。

まとめると、

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100
2					2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2		2
3					2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2		2		2	3	2
5					2		2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2
7					2		2		2	3	2		2		2	3	2		2		2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2		2	7	2	3	2		2	5	2	3	2		2		2	3	2	5	2		2	3	2		2		2	3	2	7	2		2	3	2		2	5	2	3	2		2	7	2	3	2	5	2		2	3	2

となり、最大 100 の場合でも 2 3 5 7 の倍数を埋めるだけでよかったのがわかる。

実装例

MAX = 100
invalid = Array.new(MAX.next)
(2..Integer.sqrt(MAX)).each do |e|
  unless invalid[e]
    (e * 2).step(by: e, to: MAX).each do |i|
      invalid[i] ||= e
    end
  end
end

7 は素数か？

!invalid[7]  # => true

あっている。

0 と 1 は素数か？

!invalid[0]  # => true
!invalid[1]  # => true

ちがっている。

この罠に対処するためには、

invalid[0] = true
invalid[1] = true

としてごまかしておく。

!invalid[0]  # => false
!invalid[1]  # => false

動作確認。

(0..MAX).collect { |e| e.prime? ? "o" : "." }.join     # => "..oo.o.o...o.o...o.o...o.....o.o.....o...o.o...o.....o.....o.o.....o...o.o.....o...o.....o.......o..."
(0..MAX).collect { |e| !invalid[e] ? "o" : "." }.join  # => "..oo.o.o...o.o...o.o...o.....o.o.....o...o.o...o.....o.....o.o.....o...o.o.....o...o.....o.......o..."

ヨシ！

最適化

終了判定には二通りの方法がある。

1. あらかじめ平方根を求めておく方法

(2..Integer.sqrt(MAX)).each do |e|
  # ...
end

2. 平方根をもとめず毎ループ掛け算を行う方法

(2..).each do |e|
  if e * e > MAX
    break
  end
  # ...
end

この二つはどちらが速いだろうか？ これは前者で平方根を求めた方が速い。MAX の値に関係なく速い。また MAX が大きくなるほどその差は大きくなる。

ベンチマーク

MAX = 100
require "benchmark/ips"
Benchmark.ips do |x|
  x.report("sqrt") { (2..Integer.sqrt(MAX)).each { |e| } }
  x.report("mul")  { (2..).each { |e| break if e * e > MAX } }
  x.compare!
end
# Warming up --------------------------------------
#                 sqrt   249.212k i/100ms
#                  mul   192.658k i/100ms
# Calculating -------------------------------------
#                 sqrt      2.493M (± 0.4%) i/s -     12.710M in   5.097699s
#                  mul      1.927M (± 0.4%) i/s -      9.826M in   5.097799s
#
# Comparison:
#                 sqrt:  2493283.4 i/s
#                  mul:  1927442.0 i/s - 1.29x  slower

参考

• https://www.amazon.co.jp/dp/483997750X