2015-08-06

10億までの素数を5秒で作成

改訂: 2015-09-10

• フラグを初期化していないバグがありました。
• sqrtの値(lmt)の位置の計算が含まれていませんでした。

g++ O2オプションを追加したところ、5秒以内で結果が出るようになりました。

競技プログラミング参加者は、このビットワイズ版のエラトステネスの篩を使用しているようです。 これは速いです。

kvm上のwindows vcpu 8(Xeon E3-1245)では、5秒以内で10億まで数値の素数を作り出せました。 ビルド・実行は、gccとMSYSで行いました。

説明は、下記のウェブサイトにあります。

I, ME AND MYSELF !!!: Segmented Sieve

非セグメントバージョン

sieve.cpp

#include <iostream>
#include <chrono>
#include <cmath>
#include <vector>
#include <memory>
using namespace std;
#include <stdint.h>
typedef uint32_t U32;
typedef uint64_t U64;

#define chkC(n) (pflag[n>>6]&(1<<((n>>1)&31)))
#define setC(n) (pflag[n>>6]|=(1<<((n>>1)&31)))

shared_ptr<vector<U64> > sieve(U64 n) {
  auto pprimes = make_shared<vector<U64> >(50847534L);
  vector<U64> &primes = *pprimes.get();
  U32 *pflag = (U32*)calloc(n / 64 + 1, sizeof(U32));

  U64 lmt = sqrt(n);
  U64 i, j, k;
  for (i=3; i<=lmt; i += 2) {
    if (!chkC(i)) {
      for (j=i*i, k=i<<1; j<n; j+=k) {
        setC(j);
      }
    }
  }
  primes[(j=0)++] = 2;
  for (i=3; i<n; i += 2) {
    if (!chkC(i)) {
      primes[j++] = i;
    }
  }
  primes.resize(j);
  free(pflag);
  return pprimes;
}

int main() {
  const auto startTime = chrono::system_clock::now();
  auto pprimes = sieve(1000000000L);
  const auto endTime = chrono::system_clock::now();
  const auto timeSpan = endTime - startTime;
  vector<U64>& primes = *pprimes.get();
  U32 len = primes.size();
  cout << "Prime count=" << len << ", last value=" << primes[len-1] << endl;
  cout << "Elapsed: "
    << chrono::duration_cast<chrono::milliseconds>(timeSpan)
    .count() << " ms" << endl;
  return 0;
}

ビルド

$ g++ -o sieve sieve.cpp -std=c++11 -O2

実行

$ ./sieve
Prime count=50847534, last value=999999937
Elapsed: 4849 ms