エラトステネスの篩 高速版
def _sieve(primes, min=3, size=1_0000_0000)
table = "\xFF" * (size / 8 + 1)
max = min + size * 2 - 1
tmax = size
mid = Math.sqrt(max).floor
primes.each do |i|
next if i > mid
ti = (i - min) / 2
ti += i while ti < 0
ti.step(tmax, i) do |tj|
table[tj >> 3] &= ~(1 << (tj & 7))
end
end
(0..((mid-min)/2)).each do |ti|
next unless table[ti >> 3][ti & 7].nonzero?
i = (ti << 1) + min
primes << i
ti.step(tmax, i) do |tj|
table[tj >> 3] &= ~(1 << (tj & 7))
end
end
((mid > min ? (mid - min) / 2 : 0)...tmax).each do |ti|
primes << min + ti * 2 if table[ti >> 3][ti & 7].nonzero?
end
return primes
end
def sieve(max=100, size=1_0000_0000)
max -= 1 unless max % 2 == 0
primes = []
3.step(max, size * 2) do |i|
if max > i + size * 2
primes = _sieve(primes, i, size)
else
size = (max - i) / 2 + 1
primes = _sieve(primes, i, size)
end
end
return [2]+primes
end
max = ARGV.size == 1 ? ARGV[0].to_i : 100
primes = sieve(max)
puts primes.lengthさらに高速化してみました。以下がそのテスト結果です。比較としてsuminさんの primesupto.ar.rbを載せています。
% /usr/bin/time -l ruby primesupto.ar.rb 1000000
2.12 real 2.12 user 0.00 sys
0 maximum resident set size
0 average shared memory size
0 average unshared data size
0 average unshared stack size
488 page reclaims
0 page faults
0 swaps
0 block input operations
0 block output operations
0 messages sent
0 messages received
0 signals received
0 voluntary context switches
28 involuntary context switches
% /usr/bin/time -l ruby primesupto.nr.rb 1000000
1.88 real 1.86 user 0.01 sys
0 maximum resident set size
0 average shared memory size
0 average unshared data size
0 average unshared stack size
467 page reclaims
0 page faults
0 swaps
0 block input operations
0 block output operations
0 messages sent
0 messages received
0 signals received
0 voluntary context switches
27 involuntary context switches
% /usr/bin/time -l ruby primesupto.ar.rb 10000000
24.46 real 24.40 user 0.02 sys
0 maximum resident set size
0 average shared memory size
0 average unshared data size
0 average unshared stack size
2507 page reclaims
0 page faults
0 swaps
0 block input operations
0 block output operations
0 messages sent
0 messages received
0 signals received
0 voluntary context switches
371 involuntary context switches
% /usr/bin/time -l ruby primesupto.nr.rb 10000000
19.69 real 19.67 user 0.01 sys
0 maximum resident set size
0 average shared memory size
0 average unshared data size
0 average unshared stack size
2399 page reclaims
0 page faults
0 swaps
0 block input operations
0 block output operations
0 messages sent
0 messages received
0 signals received
0 voluntary context switches
281 involuntary context switches
% /usr/bin/time -l ruby primesupto.ar.rb 100000000
427.89 real 426.75 user 0.22 sys
0 maximum resident set size
0 average shared memory size
0 average unshared data size
0 average unshared stack size
19023 page reclaims
0 page faults
0 swaps
0 block input operations
0 block output operations
0 messages sent
0 messages received
0 signals received
0 voluntary context switches
7298 involuntary context switches
% /usr/bin/time -l ruby primesupto.nr.rb 100000000
206.52 real 206.26 user 0.06 sys
0 maximum resident set size
0 average shared memory size
0 average unshared data size
0 average unshared stack size
18000 page reclaims
0 page faults
0 swaps
0 block input operations
0 block output operations
0 messages sent
0 messages received
0 signals received
0 voluntary context switches
3021 involuntary context switchesメモリ容量にあたる、page reclaims は ar も nr も同クラスになっています。これは ar で ruby の Fixnum で利用可能な 31bit のうち 8bit しか利用していないからでしょう。31bit 利用するか、nr同様 String を使うようにすれば、大きな差がつくようになると思わ
速度的には ar が予想以上に速いのに驚きます。「改良版」の1千万26秒には勝ってるし。1億になると274秒なのでだいぶ差が付きますが。
考察するに、エラトステネスの篩のホットスポットはループの最も内側である、
if (mask_bit = mask_bit_idx[idx % 2310]) != 0 then byte_idx = idx / 2310 * 60 + (mask_bit-1 << -3) mask_bit = 255 - (1 << (mask_bit & 7)) flags[byte_idx] = flags[byte_idx] & mask_bit end idx += 2 * n
table[tj >> 3] &= ~(1 << (tj & 7))
であり、この部分の速度が全体の差に大きく影響する。また、Rubyの場合、Rubyのメソッドやブロックを呼び出すと大きな差がつくが、Cのメソッドを呼ぶ分には大きな差はつかない。とすると、ar と nr は大きな差はつかないだろうと予想できます。もっとも、そう考えると今度は1億で2倍の差がつくのが疑問となります。ここで、ar の
flags = Array.new((upto + 2309) / 2310 * 60 + 60, 0xFF) # flags = "\377" * ((upto + 2309) / 2310 * 60 + 60)
のコメントアウトを逆にしてみます。すると、あら不思議・・・。
% /usr/bin/time -l ruby primesupto.ar.rb 100000000
334.33 real 333.54 user 0.25 sys
0 maximum resident set size
0 average shared memory size
0 average unshared data size
0 average unshared stack size
17114 page reclaims
0 page faults
0 swaps
0 block input operations
0 block output operations
0 messages sent
0 messages received
0 signals received
0 voluntary context switches
4913 involuntary context switches大幅に高速化されました。page reclaims も減っています。Array/int よりも String/byte の方が高速なんですねぇ。
エラトステネスの篩
なんか、エラトステネスの篩を用いていると言いつつ試行割算法を混同していたり、エラトステネスの篩を用いつつもハッシュを使ったりして速度を落としている例が多いような。エラトステネスの篩では加算しか使わない。だから速い、圧倒的に。
追記: もっとも、肝は真偽値テーブルを使うことで、一つ一つ候補を見なくて済む、篩の手法だから。割り算を排除すれば速くなるかというと、そうでもない。 (割り算がある所は高速化の余地があるというのは確かだけど)
def sieve_simple(max)
primes = []
table = Array.new(max+1).fill(true, 2);
mid = Math.sqrt(max).floor
(2..mid).each do |i|
if table[i]
primes << i
(i+i).step(max, i) do |j|
table[j] = false if table[j]
end
end
end
((mid+1)..max).each{|i|primes << i if table[i]}
primes
end多くの人のコードに存在する、n % prime == 0 が存在しないのがポイント。メモリ消費量は max+1bit なので、100万までとしても最少 数百KB となる。 (Ruby の Array だと意味ないけど、BitArray のある言語だと大きな違いになるでしょう)
ちなみに、Pentium M 1.60GHz 1GB RAM Windows XP ruby 1.9 (mingw32) で 2.298 秒でした。
参考: エラトステネスの篩 - Wikipedia
エラトステネスの篩の改良
上のソースでは探索範囲1つにつき true/false を作るので、1億まで探索しようと思うとかなりのメモリが必要となる。どうせ真偽値なのだから、BitArray を用いれば候補一つにつき 1bit で済む。また、偶数はどうせ素数ではないので、奇数のみを探索対象にすることで、テーブルサイズを半分にできる。さらに、篩範囲を大きくしていくとどうしてもマシンで確保できるメモリの限界を超えてしまうので、篩範囲を分割できるようにする。といった修正を施したバージョン。NetBSD 3.0 CeleronD 2.6GHz 512MB RAM ruby 1.9 にて 100万は数秒、1千万が1分、1億が10分でした。
参考: エラトステネスの篩い @ IDM
class BitArray
def initialize(size, bool=false)
@str = (bool ? "\xFF" : "\x00") * (size / 8 + 1)
end
def [](n)
a = n >> 3
b = n & 7
@str[a][b].nonzero?
end
def []=(n,val)
a = n >> 3
b = n & 7
if val && val != 0
@str[a] |= 1 << b if @str[a][b].zero?
else
@str[a] &= ~(1 << b)# if @str[a][b].nonzero?
end
end
def step_false(ti,tmax,i)
ti.step(tmax, i) do |tj|
a = tj >> 3
b = tj & 7
@str[a] &= ~(1 << b)# if @str[a][b].nonzero?
end
end
end
def _sieve(primes, min=3, size=1000_0000)
table = BitArray.new(size,true)
max = min + size * 2 - 1
tmax = size
mid = Math.sqrt(max).floor
primes.each do |i|
next if i > mid
ti = (i - min) / 2
ti += i while ti < 0
table.step_false(ti,tmax,i)
end
# min.step(mid, 2) do |i|
# ti = (i - min) / 2
# next unless table[ti]
(0..((mid-min)/2)).each do |ti|
next unless table[ti]
i = (ti << 1) + min
primes << i
table.step_false(ti,tmax,i)
end
((mid > min ? (mid - min) / 2 : 0)...tmax).each do |ti|
primes << min + ti * 2 if table[ti]
end
return primes
end
def sieve(max=100, size=1000_0000)
max -= 1 unless max % 2 == 0
primes = []
3.step(max, size * 2) do |i|
GC.start
if max > i + size * 2
p([:sieve, Time.now.to_f, i, size, primes.length])
primes = _sieve(primes, i, size)
else
size = (max - i) / 2 + 1
p([:sieve, Time.now.to_f, i, size, primes.length])
primes = _sieve(primes, i, size)
end
end
return [2]+primes
end追記: いろいろ高速化していたら、手元で1千万26秒、1億274秒まで短縮（’’
追記2: エラトステネスのふるいをコンパクトにする方法 に偶数以外に3や5の倍数も削る手法が。確かにこうすれば大幅にメモリ容量は減らすことができるのですが、探索する際に ti.step(tmax, i) do |tj| が使えなくなるため、大幅に速度が落ちてしまうような気もします。(偶数だけ3の倍数だけ、ならば使えるのですが)
