c++とvalarrayによるエラトステネスの篩

素数列挙プログラム

エラトステネスの篩とは、指定数以下の素数を全て見つけるためのアルゴリズムです。
詳しくはwikipediaとか見てください。

今回作るプログラムは引数の指定数以下の素数を全て列挙するものです。

いくつかの実装を作って後輩に見せびらかすアルゴリズムの工夫による速度の差を体感してもらうためにいろいろ調べていたら良いpythonのコードを見つけたので試しにvalarrayを用いてc++11互換で書きました。

ぱっと調べた限りではvalarrayとsliceを活かしたコードがなかったのでメモ感覚の投稿です。

boolの配列を作って、素数のindexだけtrueにする(素数以外をfalseにする)という仕様です。

#include <cmath>
#include <cstdint>
#include <iostream>
#include <valarray>

int main(int argc, char** argv)
{
  std::ios_base::sync_with_stdio(false);

  if (argc == 1) {
    std::cout << "Usage: " << *argv << " [number]\n" <<
                 "[number] is a positive number." << std::endl;
    return 1;
  }

  const std::uintmax_t limit {std::stoull(argv[1])}; // get search limit.
  if (limit <= 1) return 0;

  std::valarray<bool> prims(true, limit + 1); // array of [0] to [limit]
  prims[std::slice(0, 2, 1)] = false; // 0 and 1 aren't primes. (not needed)

  const std::uintmax_t search_limit {static_cast<std::uintmax_t>(std::floor(std::sqrt(static_cast<long double>(limit)))) + 1};
  for (std::uintmax_t i {2}; i < search_limit; ++i)
    if (prims[i])
      prims[std::slice(i*i, limit / i - (i - 1), i)] = false;

  for (std::uintmax_t i {2}; i <= limit; ++i)
    if (prims[i])
      std::cout << i << '\n';
}

所感

sliceが便利かつ直感的に書けて嬉しい。(実質for文が一つ減る)
環境によってはvalarrayがきっと並列処理になってくれていると信じている。

速度について

c++にvalarrayといえば気になるのが実行速度。
しかし、良い比較対象が思いつかないのでコンパイルフラグの差だけ取ります。
(vectorと取ろうかと思いましたが、for文を書くのが億劫になりました。誰かやってくd...)

以下のコードに差し替えて、篩にかける部分だけの速度を測ります。
ここからc++14互換のコードです。(auto と uniform initializationを同時に使いたいだけです)

#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <valarray>

int main(int argc, char** argv)
{
  std::ios_base::sync_with_stdio(false);

  if (argc == 1) {
    std::cout << "Usage: " << *argv << " [number]\n" <<
                 "[number] is a positive number." << std::endl;
    return 1;
  }

  const std::uintmax_t limit {std::stoull(argv[1])}; // get search limit.
  if (limit <= 1) return 0;

  const auto before_time {std::chrono::system_clock::now()};

  std::valarray<bool> prims(true, limit + 1); // array of [0] to [limit]
//  prims[std::slice(0, 2, 1)] = false; // 0 and 1 aren't primes. (not needed)

  const std::uintmax_t search_limit {static_cast<std::uintmax_t>(std::floor(std::sqrt(static_cast<long double>(limit)))) + 1};
  for (std::uintmax_t i {2}; i < search_limit; ++i)
    if (prims[i])
      prims[std::slice(i*i, limit / i - (i - 1), i)] = false;

  const auto after_time {std::chrono::system_clock::now()};

//  for (std::uintmax_t i {2}; i <= limit; ++i)
//    if (prims[i])
//      std::cout << i << '\n';

  std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(after_time - before_time).count() << std::endl;
}

実行結果はこちら。単位はミリ秒です。

$ g++ -std=c++14 -O2 main.cpp
$ ./a.out 10000000
112
$ ./a.out 10000000
111
$ ./a.out 10000000
107
$ ./a.out 100000000
1447
$ ./a.out 100000000
1455
$ ./a.out 100000000
1459

$ g++ -std=c++14 main.cpp 
$ ./a.out 100000000
2309
$ ./a.out 100000000
2477
$ ./a.out 100000000
2228

元記事のコードが10000個の要素を0.01秒で処理していたので、合わせて10000個にしようかと思いましたが 20〜65us くらいでかなりばらついたので、待ってもいいかなって思えるくらいの要素数で測りました。

要素を10倍にすると時間が約10倍強になっています。

エラトステネスの篩は計算量O(n loglogn)らしいので、オーバーヘッドがどこかにあるんだろうなーって思うとモヤモヤします。

終わりに

valarrayのsliceやgsliceなどの選択方法はとても簡素な記述ができて良いと強く思います。
もっとvalarrayを有効に使ったプログラムが増えると良いなと思います。