アトキンのふるい

説明

[0, n) の配列で，isprime[i] = true iff i は素数，な表を作る．エラトステネスのふるいよりも計算量の意味でも実用的な意味でも高速に動作する．

ただし，効率の差は N = 2^23 程度の大きな素数表に対しても高々二倍〜三倍なので，エラトステネスからこちらに切り替えたからなんとかなるという，ケースは滅多に無いと思われる．

計算量

O(n / log log n)．

実測結果

N = 16000000（≧2^23）に対して，エラトステネスのふるいとアトキンのふるいを実行したときの結果を以下に示す（30回実行したときの平均を取った）．

ソースコード

void sieve_of_atkin() {
  int n;
  for (int z = 1; z <= 5; z += 4) {
    for (int y = z; y <= sqrtN; y += 6) {
      for (int x = 1; x <= sqrtN && (n = 4*x*x+y*y) <= N; ++x)
        isprime[n] = !isprime[n];
      for (int x = y+1; x <= sqrtN && (n = 3*x*x-y*y) <= N; x += 2)
        isprime[n] = !isprime[n];
    }
  }
  for (int z = 2; z <= 4; z += 2) {
    for (int y = z; y <= sqrtN; y += 6) {
      for (int x = 1; x <= sqrtN && (n = 3*x*x+y*y) <= N; x += 2)
        isprime[n] = !isprime[n];
      for (int x = y+1; x <= sqrtN && (n = 3*x*x-y*y) <= N; x += 2)
        isprime[n] = !isprime[n];
    }
  }
  for (int y = 3; y <= sqrtN; y += 6) {
    for (int z = 1; z <= 2; ++z) {
      for (int x = z; x <= sqrtN && (n = 4*x*x+y*y) <= N; x += 3)
        isprime[n] = !isprime[n];
    }
  }
  for (int n = 5; n <= sqrtN; ++n)
    if (isprime[n])
      for (int k = n*n; k <= N; k+=n*n)
        isprime[k] = false;
  isprime[2] = isprime[3] = true;
}

分岐を消去せずに実装すると次のようになる．ただし計算量の意味でも速度は悪化する（ O(n) ）

void sieve_of_atkin() {
  int n;
  for (int x = 1; x <= sqrtN; ++x) {
    for (int y = 1; y <= sqrtN; ++y) {
      n = 4*x*x + y*y;
      if (n <= N && (n % 12 == 1 || n % 12 == 5))
        isprime[n] = !isprime[n];
      n = 3*x*x + y*y;
      if (n <= N && n % 12 == 7)
        isprime[n] = !isprime[n];
      n = 3*x*x - y*y;
      if (x > y && n <= N && n % 12 == 11)
        isprime[n] = !isprime[n];
    }
  }
  for (int n = 5; n <= sqrtN; ++n)
    if (isprime[n])
      for (int k = n*n; k <= N; k+=n*n)
        isprime[k] = false;
  isprime[2] = isprime[3] = true;
}

Verified

参考文献

• A.O.L. Atkin, D.J. Bernstein, Prime sieves using binary quadratic forms, Math. Comp. 73 (2004), 1023-1030.