アルゴリズムとデータ構造8 - 線形篩による一般乗法関数の計算

概要

競技プログラミングの問題を解く中で、線形篩を用いた一般乗法関数の計算方法を学んだので、その説明を行う。

前提知識

乗法的関数:任意の互いに素な整数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遅かった。これは定数項の違いによるものである。

タグ: 乗法的関数 線形篩 競技プログラミング 数論 アルゴリズム最適化

6月10日 19:05 投稿