区間ふるい

説明

[L, H) の配列で，table[i-L] != 0 iff i は素数，な表を作る．これは範囲を指定して作るエラトステネスのふるいであり，区間ふるいという名前がついている．素数リストだけが欲しい場合は remove--erase すればよい．

アルゴリズムは次の通り．まず高々 [0,N) の素数を普通のエラトステネスのふるいで生成する( O(N) )．次にこれらの素数を用いて [L,H) の数をふるう．区間の長さを D であらわすと計算量は O( D #(\sqrt{H} までの素数) ) である．素数定理によって H までの素数の個数は H / log(H) 程度なので，O( D \sqrt{H} / log(H) ) である．

この手法を逐次的に用いると N までの素数を空間計算量をずいぶん節約して求めることができるため，逐次ふるいと呼ばれることもある (逐次ふるいは極端には大きくない素数を列挙するために，最も現実的な方法だと考えられている)．

計算量

O( n log n )．

ソースコード

// table[i-L] == true <=> i == prime
const int SQRTN = 1<<16; // upperbound of sqrt(H)
vector<Int> segmentSieve(Int L, Int H) {
  static int p[SQRTN], lookup = 0;
  if (!lookup) {
    for (int i = 2; i < SQRTN; ++i) p[i] = i;
    for (int i = 2; i*i < SQRTN; ++i)
      if (p[i])
        for (int j = i*i; j < SQRTN; j += i)
          p[j] = 0;
    remove(p, p+SQRTN, 0);
    lookup = 1;
  }
  vector<Int> table(H - L);
  for (Int i = L; i < H; ++i) table[i - L] = i;
  for (int i = 0, j; p[i] * p[i] < H; ++i) { // O( \sqrt(H) )
    if (p[i] >= L) j = p[i] * p[i];
    else if (L % p[i] == 0) j = L;
    else                    j = L - (L % p[i]) + p[i];
    for (; j < H; j += p[i]) table[j-L] = 0;
  }
  return table;
}

逐次ふるいとしての実装．10^8 サイズの素数でも，実用的な時間で求めることができる．実際に問題に適用する場合は，これを予め計算しておいて表引きすることになる．

const int N = 100000000; // MAXPRIME 
const int M = 10000;     // SQRT(N)
const int K = 6000000;   // NUMBER OF PRIMES, CHOOSE 9/8 * N / LOG(N)
vector<Int> iterativeSieve() {
  static int p[K], table[M];
  for (int i = 2; i < M; ++i) p[i] = i;
  for (int i = 2; i*i < M; ++i)
    if (p[i])
      for (int j = i*i; j < M; j += i)
        p[j] = 0;
  p[0] = p[1] = 0;
  int num = remove(p, p+M, 0) - p;
  for (int m = M; m < N; m += M) {
    for (int x = m; x < m+M; ++x)
      table[x-m] = x;
    for (int i = 0, j; p[i]*p[i] < m+M; ++i) {
      if (p[i] >= m)          j = p[i] * p[i];
      else if (m % p[i] == 0) j = m;
      else                    j = m - (m % p[i]) + p[i];
      for (; j < m+M; j += p[i]) table[j-m] = 0;
    }
    num = remove_copy(table, table+M, p+num, 0) - p;
  }
  return vector<Int>(p, p+num);
}

Verified

• SPOJ PRIME1
• SPOJ DIVSUM2