けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

AtCoder ABC 177 E - Coprime (緑色, 500 点)

AtCoder AtCoder500点 緑色diff ABC-E エラトステネスの篩 エラトステネスの篩を用いた素因数分解の列挙 調和級数 バケット 制約:数値が10^6以下 数学(整数問題) 最大公約数 数列 最大公約数の値を固定して考える 解空間:O(N^2)通りの選択肢 ある量を固定して考える

結構教育的！！

問題概要

N 個の正の整数 A_{1}, \dots, A_{N} が与えられる。

  • これらのうちのどの 2 つも互いに素であるとき、"pairwise coprime"
  • そうではなく N 個の最大公約数が 1 であるとき、"setwise coprime
  • それ以外のとき、"not coprime"

と出力せよ。

制約

  • 2 \times N \times 10^{6}
  • 1 \times A_{i} \times 10^{6}

解法 (1)：調和級数

pairwise coprime かどうかを判定するのが本質といえる。pairwise coprime でなかったときに setwise coprime であるかどうかの判定は、「 N 個の値の GCD」を求めるだけでよい。

さて、pairwise coprime であることは、以下と同値になる。

「任意の 2 以上の整数 a に対して、 A_{1}, dots, A_{N} のうち a の倍数は高々 1 個以下である」

こうすれば解ける。あらかじめ

  • num[ v ] := A_{1}, \dots, A_{N} の中に v が何個あるか

という配列を用意しておく。そうすると、

  • 任意の a に対して
  • num[ a ], num[ 2a ], num[ 3a ], ... の総和が 1 であるかどうか

を判定すれば OK。この計算量は、 A の最大値を M として、 O(N \log M + M\log M) になる。この計算量解析には

M + \frac{M}{2} + \frac{M}{3} + \dots = O(M \log M)

という性質を用いている。

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

const int MAX = 1100000;
long long GCD(long long x, long long y) {
    if (y == 0) return x;
    else return GCD(y, x % y);
}

string solve(const vector<long long> &a) {
    vector<long long> num(MAX, 0);
    for (int i = 0; i < a.size(); ++i) num[a[i]]++;
    bool coprime = true;
    for (int a = 2; a < MAX; ++a) {
        long long sum = 0;
        for (int b = a; b < MAX; b += a) sum += num[b];
        if (sum > 1) coprime = false;
    }
    if (coprime) return "pairwise coprime";
    long long g = 0;
    for (int i = 0; i < a.size(); ++i) g = GCD(g, a[i]);
    if (g == 1) return "setwise coprime";
    else return "not coprime";
}

int main() {
    int N;
    cin >> N;
    vector<long long> a(N);
    for (int i = 0; i < N; ++i) cin >> a[i];
    cout << solve(a) << endl;   
}

解法 (2)：エラトステネスのふるいを併用して計算量改善

解法 (1) でも十分高速だが、以下の点に着目する。

「解法 (1) において、試す a は素数だけでよい」

これに着目することで、計算量を O(N \log M + M\log M) から、 O(N \log M + M \log\log M) に減らすことができる。具体的な実装としては、エラトステネスのふるいを併用することで、

  • a が素数かどうかを判定し
  • 素数 a に対してのみ、解法 (1) と同様の手続きを行う

というふうにすれば OK。

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

const int MAX = 1100000;
long long GCD(long long x, long long y) {
    if (y == 0) return x;
    else return GCD(y, x % y);
}

string solve(const vector<long long> &a) {
    vector<long long> num(MAX, 0), isprime(MAX, true);
    for (int i = 0; i < a.size(); ++i) num[a[i]]++;
    bool coprime = true;
    for (int a = 2; a < MAX; ++a) {
        if (!isprime[a]) continue;
        long long sum = 0;
        for (int b = a; b < MAX; b += a) {
            sum += num[b];
            isprime[b] = false;
        }
        if (sum > 1) coprime = false;
    }
    if (coprime) return "pairwise coprime";
    long long g = 0;
    for (int i = 0; i < a.size(); ++i) g = GCD(g, a[i]);
    if (g == 1) return "setwise coprime";
    else return "not coprime";
}

int main() {
    int N;
    cin >> N;
    vector<long long> a(N);
    for (int i = 0; i < N; ++i) cin >> a[i];
    cout << solve(a) << endl;   
}

解法 (3)：高速素因数分解

せっかくエラトステネスのふるいが出てきたので、エラトステネスのふるいを用いた高速素因数分解で解くこともできる。具体的には

  • min_prime[ v ] := v の素因子のうち最小のもの

という配列を、エラトステネスのふるいによって、以下のように求めることができる。

isprime[0] = isprime[1] = false;
min_factor[0] = 0, min_factor[1] = 1;
for (int i = 2; i <= MAX; ++i) {
     if (!isprime[i]) continue;
     min_factor[i] = i;
     for (int j = i*2; j <= MAX; j += i) {
           isprime[j] = false;
           if (min_factor[j] == -1) min_factor[j] = i;
     }
 }

これを用いると、各整数 n の素因数分解を高速に行える。具体的には 1, 2, \dots, M のすべての素因数分解が O(M \log M) で行える。

このことを用いて、以下の解法が考えられる。計算量は O(N \log M + M \log M)。

  • A_{1}, \dots, A_{N} をそれぞれ素因数分解する
  • それらが共通の素因子を持ったら "pairwise coprime" ではない

以下のコードでは、エラトステネスのふるいライブラリ (素因数分解の他、メビウス関数や約数列挙の機能も含む) を用いた。

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

long long GCD(long long x, long long y) {
    if (y == 0) return x;
    else return GCD(y, x % y);
}

// isprime[n] := is n prime?
// mebius[n] := mebius value of n
// min_factor[n] := the min prime-factor of n
struct Eratos {
    vector<int> primes;
    vector<bool> isprime;
    vector<int> mebius;
    vector<int> min_factor;

    Eratos(int MAX) : primes(),
                      isprime(MAX+1, true),
                      mebius(MAX+1, 1),
                      min_factor(MAX+1, -1) {
        isprime[0] = isprime[1] = false;
        min_factor[0] = 0, min_factor[1] = 1;
        for (int i = 2; i <= MAX; ++i) {
            if (!isprime[i]) continue;
            primes.push_back(i);
            mebius[i] = -1;
            min_factor[i] = i;
            for (int j = i*2; j <= MAX; j += i) {
                isprime[j] = false;
                if ((j / i) % i == 0) mebius[j] = 0;
                else mebius[j] = -mebius[j];
                if (min_factor[j] == -1) min_factor[j] = i;
            }
        }
    }

    // prime factorization
    vector<pair<int,int>> prime_factors(int n) {
        vector<pair<int,int> > res;
        while (n != 1) {
            int prime = min_factor[n];
            int exp = 0;
            while (min_factor[n] == prime) {
                ++exp;
                n /= prime;
            }
            res.push_back(make_pair(prime, exp));
        }
        return res;
    }

    // enumerate divisors
    vector<int> divisors(int n) {
        vector<int> res({1});
        auto pf = prime_factors(n);
        for (auto p : pf) {
            int n = (int)res.size();
            for (int i = 0; i < n; ++i) {
                int v = 1;
                for (int j = 0; j < p.second; ++j) {
                    v *= p.first;
                    res.push_back(res[i] * v);
                }
            }
        }
        return res;
    }
};

const int MAX = 1100000;
string solve(const vector<long long> &a) {
    Eratos er(MAX);
    vector<int> num(MAX, 0);
    for (auto v : a) {
        auto pf = er.prime_factors(v);
        for (auto p : pf) num[p.first]++;
    }
    bool coprime = true;
    for (int i = 0; i < MAX; ++i) if (num[i] > 1) coprime = false;
    if (coprime) return "pairwise coprime";
    long long g = 0;
    for (int i = 0; i < a.size(); ++i) g = GCD(g, a[i]);
    if (g == 1) return "setwise coprime";
    else return "not coprime";
}

int main() {
    int N;
    cin >> N;
    vector<long long> a(N);
    for (int i = 0; i < N; ++i) cin >> a[i];
    cout << solve(a) << endl;   
}