概要
競技プログラミングの問題を解く中で、線形篩を用いた一般乗法関数の計算方法を学んだので、その説明を行う。
前提知識
乗法的関数:任意の互いに素な整数p,qについてf(p)×f(q)=f(pq)を満たす関数。 完全乗法的関数:任意の整数p,qについてf(p)×f(q)=f(pq)を満たす関数。 線形篩:O(n)でn以下の素数を列挙するアルゴリズム。
適用範囲
任意の素数pに対して、f(p)およびf(p^k)を高速に計算可能な乗法的関数f。
導出過程
以下のコードは完全乗法的関数に対して正しいが、一般の場合には問題がある:
std::vector<int> primes;
std::vector<long long> func(n+1);
func[1] = 1;
for(int i=2; i<=n; ++i){
if(prime[i]){
primes.push_back(i);
func[i] = calc_prime(i); // 素数pでの値の計算
}
for(int j=0; j<primes.size(); ++j){
int64_t next = i * primes[j];
if(next > n) break;
prime[next] = false;
func[next] = func[i] * func[primes[j]]; // 乗法性を利用
if(i % primes[j] == 0) break;
}
}
この問題を解決するには、iとprimes[j]の最大公約数を除く必要がある。以下のように改良することでO(n log n)の計算量となる:
for(int j=0; j<primes.size(); ++j){
int64_t next = i * primes[j];
if(next > n) break;
prime[next] = false;
if(i % primes[j] == 0){
int a = 1, b = i;
while(b % primes[j] == 0){
a *= primes[j];
b /= primes[j];
}
if(b == 1) func[next] = calc_power(a); // 累乗の値計算
else func[next] = func[a] * func[b]; // 分解して計算
break;
}
else func[next] = func[i] * func[primes[j]];
}
さらに最適化するために、最小素因数の累乗を記録するlow配列を導入する:
std::vector<int> low(n+1);
func[1] = 1; low[1] = 1;
for(int i=2; i<=n; ++i){
if(prime[i]){
primes.push_back(i);
low[i] = i;
func[i] = calc_prime(i);
}
for(int j=0; j<primes.size(); ++j){
int64_t next = i * primes[j];
if(next > n) break;
prime[next] = false;
if(i % primes[j] == 0){
low[next] = low[i] * primes[j];
int a = low[next], b = next / a;
if(b == 1) func[next] = calc_power(a);
else func[next] = func[a] * func[b];
break;
} else {
low[next] = primes[j];
func[next] = func[i] * func[primes[j]];
}
}
}
興味深い事実
LOJ #124の問題において、理論上はO(n)のアルゴリズムがO(n log n)より速いはずだが、実際の実行時間では1039ms遅かった。これは定数項の違いによるものである。