エラトステネスの篩の省メモリ化

はじめに

$n$ 未満の素数を列挙する方法として、エラトステネスの篩があります．
エラトステネスの篩は計算に $O(n)$ ビット領域必要とします．
ここでは、区間篩を用いて $O\left(\sqrt{n}\right)$ ビット領域に減らせることを示します．

1. エラトステネスの篩

エラトステネスの篩は $n$ 未満のすべての素数を列挙します．次の手順で計算します．

1. $2$ 以上 $n$ 未満の数を表に並べます．

2. $i=2,3,\ldots,\left\lceil\sqrt{n}\right\rceil-1$ について、以下の操作を繰り返します．

   • 表に $i$ が存在するとき、表から $i^2$ 以上の $i$ の倍数を取り除く．

3. 最終的に残った数が $n$ 未満の素数となります．

実装は次のようになります．

sieve.cpp

vector<int> sieve(int n) {
    vector<bool> is_prime(n, true);
    is_prime[0] = is_prime[1] = false;
    int sqrt_n = (int)ceil(sqrt(n));
    for (int i = 2; i < sqrt_n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j < n; j += i) {
            is_prime[j] = false;
        }
    }
    vector<int> primes;
    for (int i = 2; i < n; ++i) {
        if (is_prime[i]) primes.push_back(i);
    }
    return primes;
}

上のコードの計算量は $O(n\ln{\ln{n}})$ 時間、$O(n+\pi(n)\lg{n})$ ビット領域です．ただし、$\pi(n)$ は $n$ 以下の素数の個数です．素数フラグ is_prime が $O(n)$ ビット領域、素数リスト primes が $O(\pi(n)\lg{n})$ ビット領域です．

2. 区間篩

区間篩は $a$ 以上 $b$ 未満のすべての素数を列挙します．手順は次の通りです．

1. エラトステネスの篩を用いて、$1$ 以上 $\sqrt{b}$ 未満の素数列 $p_1,p_2,\ldots,p_m$ を計算します．

2. $a$ 以上 $b$ 未満の数を表に並べます．

3. $1\le i\le m$ について、$a$ 以上 $b$ 未満の $p_i$ の倍数を表から取り除きます．

4. 最終的に残った数が $a$ 以上 $b$ 未満の素数となります．

実装は次のようになります．

segment-sieve.cpp

vector<int> segment_sieve(int a, int b) {
    int sqrt_b = (int)ceil(sqrt(b));
    vector<bool> is_prime_small(sqrt_b, true);
    is_prime_small[0] = is_prime_small[1] = false;
    for (int i = 2; i < sqrt_b; ++i) {
        if (!is_prime_small[i]) continue;
        for (int j = i * i; j < sqrt_b; j += i) {
            is_prime_small[j] = false;
        }
    }

    vector<bool> is_prime(b - a, true);
    for (int i = 2; i < sqrt_b; ++i) {
        if (!is_prime_small[i]) continue;
        int k = max(i, (a + i - 1) / i);
        for (int j = k * i; j < b; j += i) {
            is_prime[j - a] = false;
        }
    }

    vector<int> primes;
    for (int i = 0; i < b - a; ++i) {
        if (is_prime[i] && i + a >= 2) {
            primes.push_back(i + a);
        }
    }
    return primes;
}

上のコードの計算量は $O((N+\sqrt{b})\ln{\ln{b}})$ 時間、$O(N+\sqrt{b}+\pi(N)\lg{N})$ ビット領域です．ただし、$N=b-a$ です．

3. 本題

エラトステネスの篩を用いて $n$ 未満の素数を $O(n)$ ビット領域で計算することができました．区間篩を用いることで $O\left(\sqrt{n}\right)$ ビット領域に減らすことができます．手順は次の通りです．

1. 区間 $[0,n)$ を長さ $\left\lceil\sqrt{n}\right\rceil$ の小さい区間に分割します．
   分割した各区間を $[a_1,b_1),[a_2,b_2),\ldots,[a_m,b_m)$ とします．
   このとき、$a_1=0,b_s-a_s=\left\lceil\sqrt{n}\right\rceil;(1\le s\le m)$ を満たします．

2. $s=1,2,\ldots,m$ について、$[a_s,b_s)$ に対して区間篩を行います．

実装は次のようになります．

efficient-space-sieve.cpp

vector<int> efficient_space_sieve(int n) {
    int sqrt_n = (int)ceil(sqrt(n));
    vector<bool> is_prime_small(sqrt_n, true);
    is_prime_small[0] = is_prime_small[1] = false;
    for (int i = 2; i * i < sqrt_n; ++i) {
        if (!is_prime_small[i]) continue;
        for (int j = i * i; j < sqrt_n; j += i) {
            is_prime_small[j] = false;
        }
    }

    vector<int> primes;
    int m = (n + sqrt_n - 1) / sqrt_n;
    for (int s = 1; s <= m; ++s) {
        int a = (s - 1) * sqrt_n, b = s * sqrt_n;
        vector<bool> is_prime(b - a, true);
        for (int i = 2; i < sqrt_n; ++i) {
            if (!is_prime_small[i]) continue;
            int k = max(i, (a + i - 1) / i);
            for (int j = k * i; j < b; j += i) {
                is_prime[j - a] = false;
            }
        }
        for (int i = 0; i < b - a; i++) {
            if (is_prime[i] && i + a >= 2 && i + a <= n) {
                primes.push_back(i + a);
            }
        }
    }
    return primes;
}

上のコードの計算量は $O(n\ln{\ln{n}})$ 時間、$O\left(\sqrt{n}+\pi(n)\lg{n}\right)$ ビット領域です．

4. 計測

計測用に下のコードを使用します．比較のため、各関数は素数リストではなく、素数の個数を計算しています．sieve1.cpp, sieve2.cpp の領域計算量はそれぞれ $O(n),O\left(\sqrt{n}\right)$ ビット領域です．

sieve1.cpp

#include <bits/stdc++.h>
using namespace std;

int sieve(int n) {
    vector<bool> is_prime(n, true);
    is_prime[0] = is_prime[1] = false;
    int sqrt_n = (int)ceil(sqrt(n));
    for (int i = 2; i < sqrt_n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j < n; j += i) {
            is_prime[j] = false;
        }
    }
    int cnt = 0;
    for (int i = 2; i < n; ++i) {
        if (is_prime[i]) ++cnt;
    }
    return cnt;
}

int main() {
    int n;
    cin >> n;
    cout << sieve(n) << endl;
}

sieve2.cpp

#include <bits/stdc++.h>
using namespace std;

int efficient_space_sieve(int n) {
    int sqrt_n = (int)ceil(sqrt(n));
    vector<bool> is_prime_small(sqrt_n, true);
    is_prime_small[0] = is_prime_small[1] = false;
    for (int i = 2; i * i < sqrt_n; ++i) {
        if (!is_prime_small[i]) continue;
        for (int j = i * i; j < sqrt_n; j += i) {
            is_prime_small[j] = false;
        }
    }

    int cnt = 0;
    int m = (n + sqrt_n - 1) / sqrt_n;
    for (int s = 1; s <= m; ++s) {
        int a = (s - 1) * sqrt_n, b = s * sqrt_n;
        vector<bool> is_prime(b - a, true);
        for (int i = 2; i < sqrt_n; ++i) {
            if (!is_prime_small[i]) continue;
            int k = max(i, (a + i - 1) / i);
            for (int j = k * i; j < b; j += i) {
                is_prime[j - a] = false;
            }
        }
        for (int i = 0; i < b - a; i++) {
            if (is_prime[i] && i + a >= 2 && i + a <= n) {
              	++cnt;
            }
        }
    }
    return cnt;
}

int main() {
    int n;
    cin >> n;
    cout << efficient_space_sieve(n) << endl;
}

$N=10^6,10^7,10^8,5\times10^8,10^9$ における計測結果は次の通りです．

$N$	sieve1.cpp	sieve2.cpp (省メモリ)
$10^6$	11.6 ms / 560.8 KB	14.4 ms / 430.4 KB
$10^7$	51.2 ms / 1588.0 KB	80.8 ms / 449.6 KB
$10^8$	460.4 ms / 12568.8 KB	672.2 ms / 446.4 KB
$5\times10^8$	3727.6 ms / 61524.8 KB	3260.8 ms / 452.0 KB
$10^9$	8438.2 ms / 122684.8 KB	6451.6 ms / 451.2 KB

全体的に sieve2.cpp の方がメモリ使用量が少ないことが分かります．
また、$N$ が $10^8$ 以下のときは sieve1.cpp の方が速いですが、$N$ が $5\cdot 10^8$ 以上の時は sieve2.cpp の方が速くなっています．

参考