C++によるエラトステネスの篩の実装

こんにちは、Sayahamittです。

2月中旬は本州太平洋側が二度の記録的な大雪に見舞われ、我が家も慣れない雪かきに追われました。

雪かきをしてしまった後、やることもないので久々に素数に遊ばれようかと思った次第です。(とても素数で遊んだ、と言えるレベルじゃない(´・ω・｀)

さて、素数を求めるアルゴリズムは幾通りもありますが、今回は素数探索アルゴリズムの代表格であるエラトステネスの篩を素直に、バカ正直に実装してみたので覚書程度に記事にします。

C++で、と銘打ってありますが楽をするためにC++を使っただけです。本当はC言語で書くべきである程度のコードしか書いていませんがご容赦ください。

エラトステネスの篩

エラトステネスの篩は素数を探索したい範囲を0-nとして予め決めておき、その中で既知の素数を因数に持つ数をふるい落とすというものです。ここではソースコードを書くという視点から注意すべき(私が注意した)ポイントを中心にエラトステネスの篩をザッと説明しようと思います。

手順

Nまでの素数リストを得るには、0〜Nまでの整数のリストを作る。
整数リストから0と1を除く。
最初の素数である2を素数リストに移し、整数リストから2の倍数、つまり「2を因数に持つ数」を除く(2で篩う)
整数リストの先頭(2は素数リストへ移したので現在は3)を素数リストへ移し、その数の倍数を整数リストから除く。(3で篩う)

以降、手順4を整数リストの先頭の数の平方がNより大きくなるまで繰り返す。

終了条件 : 整数リストの先頭の数の平方がNより大きい。

このアルゴリズムのミソは、

素数の倍数のみで整数リストを篩えばよいこと。

ある数で整数リストを篩ったあと、整数リストの次の数は必ず素数であること。

Nまでの素数を求めたい時N^(1/2)までの数でリストを篩えばよいので計算の手間が比較的少なくて済むこと

です。

実装

これをbool型の配列一つを用いて実装します。

整数リストを配列のインデックス(添字)とし、ある要素についてインデックスが素数であればtrue、そうでなければfalseとします。

エラトステネスの篩では「整数リスト」と「素数リスト」が必要になるわけですが、実際に計算を行う時に配列を2つ作ると素数をコピーする手間が増えますし、そもそもメモリの無駄です。

そこで、N以下の素数リストを求めたい時にはN個の要素を持つ配列を作り、そのインデックス(要素番号)が0からN-1になる事を利用して、ある数mが素数かどうかをm番の要素に記録することにします。

すると、例えば3までで篩った後には2と3の倍数の要素にはその番号が素数でない旨が”false”として記録されています。

この状態で整数リストの先頭を得るには、”false”と記されている要素番号を読み飛ばして配列を4番目の要素から一要素ずつチェックして、最初に出てきた”true”の要素番号を取れば良いのです。

これをfor文で実装すると以下のようになります。

void SieveofEratosthenes(long long askmax){
  //引数askmaxで指定された大きさの配列を動的確保
  //(askmaxまでの素数リストを求める)
  bool* seq = new bool[askmax];

  //配列をtrueで初期化
  fill(seq,seq+askmax,true);

  //0と1は例外的に非素数として処理、予めfalseとしておく
  seq[0]=seq[1]=false;

  //ここからエラトステネスの篩
  //numを整数リストの先頭の数とする(リストを篩う数)
  for(long long num=2;num*num<askmax;num++){

    //numが整数リストの先頭なら、つまり素数ならリストを篩う
    if(seq[num]==true){
      //numの倍数をfalseにすれば良いのでnum個目毎の要素を非素数としてfalseにする
      for(long long cursor=num*num;cursor<askmax;cursor += num){
	seq[cursor]=false;
      }
    }
  }

  delete seq;
}

このコードでは一段目のforで配列seqを一要素ずつチェックし、

tureの要素が出てきたらそのインデックスの数で配列(整数リスト)を篩うようにしています。

エラトステネスの篩の二段目のfor文(20行目)について少し説明します。

倍数を求めるのに割った余りが0であるか調べなくて良いのか?

ある数mがnumの倍数かどうかを調べようと思うと、つい

if( m%num == 0 )

と書きたくなってしまうのですが、配列の中でnumの倍数である要素番号というのはnum個ごとに一つあるわけで、for文を用いてnum個ごとに配列の要素を参照することは

for(long long cursor=foo; bar; cursor += num)

と書けば容易にできるわけです。

「篩う数」の2乗の数から整数リストを篩い始めているのはなぜか

二段目のfor(20行目)で配列seqへのアクセス子として使っているcursorをnum*numで初期化しているのは何故かと言うと
num で篩おうとしている時には、 num 以下の素数で既に整数リストが篩われているからです。

m<num とした時に(m は num 以下の任意の素数) m*num は num の倍数でもありますが、同時に m の倍数であり、 num より小さい素数 m で既に篩われているのです。

全ての合成数は複数の素数の乗算として表されるので、以上のことは m<num の場合 num*num 未満の全ての合成数に当てはまり、結果として num で整数リストを篩う時には num*num の数から始めれば良いことになります。

上記のエラトステネスの篩を動かす為の完動コードも一応書いておきます。

#include <iostream>
#include <time.h>

void SieveofEratosthenes(long long);
long long showprim(bool*,long long,bool);

int main(){
  long long input=0;
  //整数リストの大きさを入力させる
  std::cout<<"求める素数の最大値を入力してください >"<<std::flush;   std::cin>>input;
  SieveofEratosthenes(input);
  return 0;
}

//整数リスト配列のうち素数だけを表示、最大の素数を返す関数
//引数modeに1を与えると与えられた配列中の最大の素数を返すだけ
//modeに0を与えると配列中の全ての素数を列挙表示する。
long long showprime(bool* seq,long long askmax,bool mode){
  long long lpn=0;
  int i=askmax-1;
  int j=0;

  if(mode==0){
    for(long long i=0;i<askmax;i++){
      if(seq[i]==1){
	std::cout<<i<<", "<<std::flush;
      }
    }
    std::cout<<"over"<<std::endl;
  }

  //配列の先頭から最大素数を探すと時間が掛かるので
  //i=askmax-1としておいて配列の末尾から素数を探して最初の素数を返す。
  do{
    if(seq[i]==true){
      lpn=i;
    }
    i--;
  }while(lpn==0);

  return lpn;
}

void SieveofEratosthenes(long long askmax){
  //時間計測用変数の宣言
  time_t start,end;

  bool* seq = new bool[askmax];
  std::fill(seq,seq+askmax,true);
  seq[0]=seq[1]=false;

  time(&start);
  //---エラトステネスの篩---
  for(long long num=2;num*num<askmax;num++){
    if(seq[num]==true){
      for(long long cursor=num*num;cursor<askmax;cursor += num){
	seq[cursor]=false;
      }
    }
  }
  //-----------------------
  time(&end);
  std::cout<<"runtime : "<<end-start<<"(sec)"<<std::endl
  <<"Largest prime number : "<<showprime(seq,askmax,1)<<std::endl;

  delete seq;
}

このコードでは1億までの素数リストを求めるのに12秒、10億までで158秒かかります。
(実験環境:core i7 3820 3.6Ghz, RAM DDR3-1333 4channel 16Gbyte)

他のサイトやブログで紹介されているものと比べると遅いので改良の余地がありそうです。

また、ユーザーからの入力に応じて愚直に配列を動的確保しているので、実行しているマシンで確保できるメモリ容量を超える大きさの素数リストを求めようとするとコアダンプエラーでOSに落とされます。

現状ではクズコードですが、配列の分割やマルチスレッド化も視野に入れつつ気長に弄って行こうと思います。

エラトステネスの篩

実装

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル