線形ふるい法 (Linear Sieve)
線形ふるい法は、ある範囲内の素数を効率的に求めるアルゴリズムです。このアルゴリズムの時間計算量はO(n)であり、各合成数が一度だけふるい落とされることを保証します。本稿では、線形ふるい法を用いて素数、オイラー関数、メビウス関数を計算する方法を解説します。
1. 素数の線形ふるい
素数の線形ふるいは、指定された上限Nまでのすべての素数をO(N)の時間で求めるアルゴリズムです。その核心は、各合成数がその最小の素因数によってのみふるい落とされるという点にあります。これにより、重複してふるい落とす無駄がなくなり、線形時間で処理が完了します。
#include <iostream>
#include <vector>
#define MAX_N 10000050
/*
線形ふるい法
O(n)でn以下の素数を列挙する
*/
int main() {
std::vector<bool> isComposite(MAX_N, false);
std::vector<int> primes;
int n;
std::cin >> n;
for (int i = 2; i <= n; ++i) {
if (!isComposite[i]) {
primes.push_back(i);
}
for (int p : primes) {
if (i * p > n) {
break;
}
isComposite[i * p] = true;
if (i % p == 0) {
break;
}
}
}
// primes配列に素数が格納されている
return 0;
}
if (i % p == 0) break; という条件がこのアルゴリズムの鍵です。この条件により、合成数i*pは、iの最小の素因数pによってのみふるい落とされます。もしこの条件がない場合、合成数は複数の素因数の組み合わせによって複数回ふるい落とされてしまい、計算量がO(n log log n)に悪化します。
2. オイラー関数の線形ふるい
オイラー関数φ(n)は、1以上n未満の整数のうち、nと互いに素なものの個数を表します。線形ふるい法は、素数ふるいと同様の構造を利用して、オイラー関数の値も効率的に計算できます。
#include <iostream>
#include <vector>
#define MAX_N 10000050
/*
線形ふるい法
O(n)でn以下の素数とオイラー関数の値を求める
*/
int main() {
std::vector<bool> isComposite(MAX_N, false);
std::vector<int> primes;
std::vector<int> eulerPhi(MAX_N);
int n;
std::cin >> n;
eulerPhi[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!isComposite[i]) {
primes.push_back(i);
eulerPhi[i] = i - 1;
}
for (int p : primes) {
if (i * p > n) {
break;
}
isComposite[i * p] = true;
if (i % p == 0) {
eulerPhi[i * p] = eulerPhi[i] * p;
break;
} else {
eulerPhi[i * p] = eulerPhi[i] * (p - 1);
}
}
}
return 0;
}
このアルゴリズムは、素数ふるいのループ内でオイラー関数の値を更新します。もしiが素数pで割り切れる場合（iの最小の素因数がpである場合）、φ(i*p) = φ(i) * p となります。そうでない場合、iとpは互いに素であるため、φ(i*p) = φ(i) * φ(p) = φ(i) * (p - 1) となります。
3. メビウス関数の線形ふるい
メビウス関数μ(n)は、数論における重要な積性関数の一つです。その定義は以下の通りです。
- μ(1) = 1
- もしnが平方因子を持つなら、μ(n) = 0
- もしnがk個の異なる素因数を持つ素数の積であるなら、μ(n) = (-1)^k
線形ふるい法を用いて、メビウス関数の値も同様に計算できます。
#include <iostream>
#include <vector>
#define MAX_N 10000050
/*
線形ふるい法
O(n)でn以下の素数とメビウス関数の値を求める
*/
int main() {
std::vector<bool> isComposite(MAX_N, false);
std::vector<int> primes;
std::vector<int> mobius(MAX_N);
int n;
std::cin >> n;
mobius[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!isComposite[i]) {
primes.push_back(i);
mobius[i] = -1;
}
for (int p : primes) {
if (i * p > n) {
break;
}
isComposite[i * p] = true;
if (i % p == 0) {
mobius[i * p] = 0;
break;
} else {
mobius[i * p] = -mobius[i];
}
}
}
return 0;
}
このアルゴリズムも同様の構造を持ちます。もしiが素数pで割り切れる場合、i*pにはp^2という平方因子が含まれるため、μ(i*p) = 0 となります。そうでない場合、i*pの素因数の個数はiの個数と1増えるため、符号が反転し、μ(i*p) = -μ(i) となります。