エラトステネスの篩い

最終更新日:2003.05.31

IDM | 素因数分解 | エラトステネスの篩い

目次

  1. 紹介
  2. アルゴリズムの基本
  3. Rubyによる実装
  4. 改良1:偶数を省く
  5. 改良2:配列の分割
  6. 実装
  7. 関連資料
  8. 変更履歴

紹介

最も古典的な素数表作成アルゴリズムとして、エラトステネスの篩い(Eratosthenes's sieve)があります。

今ではより効率的な代替アルゴリズムもありますが、次に挙げる理由から、ここではエラトステネスの篩いを説明します。

一番上に戻る

アルゴリズムの基本

  1. 求める素数の上限値をMとします。初めに、2からMまでの整数を全て書きます。

    計算機に於いては大きさ(M-1)の配列を確保して0で埋めれば十分でしょう。

  2. 次にその整数リストの中から、2以外の2の倍数全てに対し、「合成数」マークを付けます。

    手計算を行うときはその数字の上に斜線を引いたりしますが、計算機に於いては配列の対応する要素に1を代入すればよいでしょう。

  3. 次に、2より大きく「合成数」マークの付いていない最小の整数を探します。当然の事ながら、3が見つかります。

  4. 3自身以外の3の倍数全てに対し、「合成数」マークを付けます。

  5. 以下同様にして、順次「合成数」マークの付いていない最小の整数pを探して、その倍数に「合成数」マークを付けていきます。

pがM0=√Mを超えるまでこれを繰り返せば、マークが付かずに残った整数は素数です。

合成数マークを付けることを篩いに喩え、マークが付かずに残った素数を「篩の目をくぐった」と考えるので、「エラトステネスの篩い」と言います。

証明

2からMまでの整数に対して、「マークされていないこと」と「素数であること」が同値であることを証明する。

  1. 2は素数であり、上記アルゴリズムでは決してマークされない。従って、2に対しては上記命題が成立する。

  2. n未満(n≦M)の任意の整数に対して上記命題が成立すると仮定する。

    nがマークされたならば、nはn以外のある数pの倍数である。上記アルゴリズムではp=1では有り得ないから、このときnは合成数である。

    nが合成数とすると、nはn未満の素因数pをもつ。仮定よりpはマークされていないので、上記アルゴリズムに従ってpの倍数はマークされる。従ってnはマークされる。

    従って、nに対しても上記命題が成立する。

  3. 数学的帰納法により、2からMまでの任意の整数に対して「マークされていないこと」と「素数であること」は同値である。

一番上に戻る

Rubyによる実装


#!/usr/local/bin/ruby

m = ARGV.shift.to_i
WorkArea = Array.new(M+1, false)   # 配列を確保してfalseで埋める

m_0 = Math.sqrt(m).to_i
p=2
while p <= m_0 do
  n=p
  while (n+=p) <= m do
    WorkArea[n] = true         # pの倍数に合成数マークをつける
  end

  begin
    p += 1                     # 次の最小未マーク整数を探す
  end while WorkArea[p]
end

2.upto(m) do |n|
  print "#{n}\n" unless WorkArea[n] # マークされていなければ出力
end

一番上に戻る

改良1:偶数を省く

2が素数であること、それ以外の偶数が合成数であることはあまりにも明らかです。しかし、上記実装ではそのことをわざわざ計算して求めています。

範囲内の整数全てを素数候補とするのではなく、奇数のみを候補にすれば計算量と必要記憶容量は半減します。これは素朴試し割り法のときと同じです。


#!/usr/local/bin/ruby

m = ARGV.shift.to_i
workarea_size = (m-1)/2+1

# WorkArea[i]は整数2*i+3に対応する
WorkArea = Array.new(workarea_size, false)  


i_max = (Math.sqrt(m).to_i - 1)/2    # ループの上限値
i=0
while i<=i_max do
  p=2*i+3
  j=i
  while (j+=p)<=workarea_size do
    WorkArea[j] = true          # 倍数に対応する要素をマークする
  end
  begin
    i+=1                        # 次の最小未マーク整数を探す
  end while WorkArea[i]
end

print "2\n"
WorkArea.each_index do |i|
  print "#{2*i+3}\n" unless WorkArea[i]
end

同様にして3,5,...の倍数を予め除いておくことも可能ではあるでしょうが、アルゴリズムが複雑になるので今は採用しません。

一番上に戻る

改良2:配列の分割

上記実装では、M以下の奇数に対応する要素を、一挙に１つの配列として確保しています。しかし、これでは2の32乗以下の素数を求めるだけでも膨大なメモリーを必要とします。C言語でchar型配列を用いたとしても2GBのメモリーを必要とすることになり、（仮想メモリーを使用するにせよ）2002年現在一般的な性能のパソコンでは実行が難しくなってきます。

実際、改良1を施したC言語による実装を手元の環境で走らせてみたところ、上限として2の29乗を与えたところで落ちました。

素数定理より、上限を大きくするほど整数に占める素数の割合が少なくなっていくことを考えると、上限に比例するメモリーを要求することはなんとも虚しい気がします。

さて、この点を何とかしましょう。

合成数Nは√N以下の素因数を持っているため、Nに「合成数」マークを与えるためには√N以下の素数表が得られていれば十分です。そこで、Mまでの素数表を作成するために、まず√Mまでの範囲に対して上記と同様の篩を行って素数表を作ます。そして、(√M)+1からMまでの範囲に対してはその小素数表に載っている素数の倍数をマークすることにすれば十分です。

このことにより、√M+1からMまでの部分については、何回かに分けて「篩う」ことが容易になります。


#!/usr/local/bin/ruby

m = ARGV.shift.to_i

print "2\n" if m>=2

# ワークエリアのサイズを決定
#  mに対してあまり細かすぎたり、絶対的に小さすぎたりすると
#  ループ回数が増え、効率が悪い
#  かと言って、あまり大きすぎるとメモリアクセスの関係から
#  効率が悪い。
N_SPLIT = 10
N_WS_MAX = 256*1024  
N_WS_MIN =  10*1024
workarea_size = (m/N_SPLIT).to_i
if workarea_size < N_WS_MIN then
  workarea_size = N_WS_MIN
elsif workarea_size > N_WS_MAX then
  workarea_size = N_WS_MAX
end

# WorkAreaを確保
WorkArea = Array.new(workarea_size)


# m_0 = [sqrt(m)]以下の素数表を得る
# WorkArea[i] は整数2*i+3に対応する

SmallPrimeTable = Array.new

# 小素数表の上限
m_0 = Math.sqrt(m).to_i

# i_max : 2*i+3 <= sqrt(m_0)なる最大のi
i_max = ((Math.sqrt(m_0).to_i-3)/2).to_i
# j_max : 2*j+3 <= m_0 なる最大のj
j_max = ((m_0 - 3)/2).to_i


# このあたりは改良1のときと同じ
i = 0
while i <= i_max do
  p = 2*i+3
  j = i
  while (j+=p) <= j_max do
    WorkArea[j] = true
  end

  begin
    i += 1
  end while WorkArea[i]
end

SmallPrimeTable = Array.new

0.upto(j_max) do |i|
  unless WorkArea[i] then
    SmallPrimeTable.push 2*i+3  # 未マークならば、素数表に追加。
    print "#{2*i+3}\n"
  end
end


# (2*i+base)%p==0 なる最小の非負整数iを与える関数
def least_index(base, p)
  i=-base%p
  i+=p if i<0
  if i%2==0 then
    return i/2
  else
    return (i+p)/2
  end
end


# 残りのm_0...m を篩にかける
# WorkArea[i] は整数2*i+baseに対応する

base = (2*j_max+3)+2  # m_0より大きい最小の奇数

while base < m do
  # WorkAreaを初期化
  WorkArea.fill(nil)

  # n_max = min(base+2*workarea, m) に対して、
  # 2*i+base <= n_max  なる最大のiをi_maxとする。
  #（ループカウンタiの上限)。
  i_max = if 2*(workarea_size-1)+base > m then
	    (m-base)/2.to_i
	  else
	    workarea_size-1
	  end

  # 小素数表の各素数pに対して、
  SmallPrimeTable.each do |p|
    # p|(2*i+base)なる最小のiを与える
    i = least_index(base, p)
    while i <= i_max do
      WorkArea[i] = true        
      i += p
    end
  end
  
  0.upto(i_max) do |i|
    # 未マークならば、出力する。
    print "#{2*i+base}\n" unless WorkArea[i]
  end

  # 次のbaseを算出。
  base += 2*workarea_size
end

一番上に戻る

実装

実装一覧
言語所在サイズ(bytes)備考
rubyimpl/eratosthenes/esieve.rb3,029
Cimpl/eratosthenes/esieve2.c4,243整数はunsigned int型を使っている。 MAX_UINT以下の範囲しか篩えない。
JAVAimpl/eratosthenes/esieve.java4,041
上記全てのアーカイブimpl/eratosthenes/esieve.tar.gz4,613
JavaScriptimpl/eratosthenes/esieve.js2345

一番上に戻る

関連資料

一番上に戻る

変更履歴

2002.11.01

関連資料の項を追加。

3,5,7の倍数を予め除いてもあまり効率が上がらなかったのは、アルゴリズムとしての効率の問題ではなく、単に私の実装が悪かっただけの模様。そのあたりの記述を訂正した。そのうち、ちゃんと再実装したいと思う。

2002.09.22

なんか、「1〜1000までの素数」というキーワードでこのページへやって来る人が多いです。

この人たちの中には、C,JavaのコンパイラもRubyインタプリタも手元にない人もいるであろうと考えて、JavaScript版を作ってみました。

2002.08.26

文中のRubyによる実装を（寝ぼけて）いじくり回しているうちに、整合性がとれなくなっていたらしい。バグを修正。

2002.07.29

証明の記述が変だったのを修正。

2002.07.24

C言語実装esieve2.cにて、32bit以下の素数表作成に成功しました。次の環境で所要時間は約1時間、出力は813,120,884bytesでした。

機種Dell Dimension 4100
CPUPentium III 800MHz
RAMSDRAM 256MB
OSWindows Me + Cygwin
使用コンパイラgcc 2.95.3-5 (with -O3 option)

但し、実際に用いたのは、出力を標準出力でなくバイナリファイルに切り替えたesieve4.cです（esieve2.cそのままの場合、画面出力に時間を食うのでもっと遅いと思われます）。カレントディレクトリのprime.datに出力します。

2002.07.24

小素数表作成のバグを修正。sqrt(sqrt(m))以下の素数しか列挙していなかった。

一番上に戻る

IDM | 素因数分解 | エラトステネスの篩い

© 2002-2003 ゆうき