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

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

AtCoder AGC 003 D - Anticube (1100 点)

とくらちゃん・シンヤカトー・てんぷらたんたちと気合いで解いた。

問題へのリンク

問題概要

 N 個の正の整数  s_1, s_2, \dots, s_N が与えられる。この中から最大個数の要素を取り出して

  • どの 2 つをとってもその積が立方数とはならない

という条件を満たすようにせよ。

制約

  •  1 \le N \le 10^{5}
  •  1 \le s_i \le 10^{10}

解法

パッと見て一応「最大独立集合問題」に属する問題ではある。

素因数分解的なことを考えたくなるので、とりあえず s たちが素数べきの場合を考えてみる。そうすると、s たちは素数べきの指数の mod. 3 のみに着目すればよくて

  • 0 のやつは複数あっても 1 個しか採用できないし、他と組み合わさって 0 になることはないから、とりあえず 1 個除外して考えてしまってよい
  • 1 のやつと 2 のやつとがペアになる (なので、1 のやつと 2 のやつのうち、個数が多い方が答え)

という感じになっていることがわかる。

そして (最終的にこの考察はあまり意味がなかったが)、s たちを一般の整数にしても、予め立方数を取り除いておくと、「このペアは立方数になる」という関係でグラフを作ると二部グラフになることがわかった。理由は奇数長のサイクルを構成しえないから。一瞬、二部グラフの最大独立集合問題じゃーーーーん!!!とぬか喜びをした。次の瞬間、グラフが dense になりうることに気がついて、ダメじゃーーーんとなった。

次にふと制約に目をやるとそもそも

  • 全部の s を素因数分解するのが間に合わない!!!!!!!!!!!

ということに気がつく!!!!!! え!?となり、でも「素数の指数」のみが重要だから、105 までの素数についてのみ割れば OK と思い、でもそれでも TLE する計算になった。

こうなったら、 10^{3.3} \le 2500 で、2500 以下の素数までで良いようになる解法があるに違いない、そういう制約に違いないと思って考察を頑張った。

まず、素因数分解したときの各指数 mod. 3 のみが重要なので、mod. 3 したくなる。そしてそれは確かに 2500 以下の素数について指数を mod. 3 すれば実現できる。2500 より大きい素数については、指数が 2 より大きくはならないからだ。というわけで、すべての s に対して、指数が mod. 3 になるように正規化しておくことにする。

これに気がつくと、ある重要な性質が浮かび上がった。

  • ab と ac が立方数ならば、b = c

が成り立つということだった。つまり、a と掛け算して立方数になれる相手は一意に定まるのだ!!!!!!!!!

ここまで来ると我々の考えるべき道筋はただ 1 つ。


いかにして、a が与えられて、ab が立方数になるような b を高速に求められるか???


制限として、a を完全な形で素因数分解するのでは間に合わない。しかし、a を 2500 以下の素数について試し割りすることはできる。

 a = p_1^{x_1} ... p_K^{x_K} (残り) ( x_i = 1 or 2)

という風にしたときに、とりあえず

 b = p_1^{3 - x_1} ... p_K^{3 - x_K} (残り)

という形になることはわかる。(残り) のところを頑張る。実は  a についての (残り) のあり方は以下の 4 パターンしかない (素数が 3 個以上になると  10^{10} を超えるため)

1 個目のパターンは自明。3 個目のパターンになるかどうかは判定できる。そして

  • 3 個目、相方は P
  • 2 個目と 4 個目は、どちらであっても相方は二乗したもの

であることがわかる。これで全体を通してなんとか解けそうとなった。

#include <iostream>
#include <vector>
#include <map>
using namespace std;

const int MAX = 2500;      // 素数判定する最大の数
bool IsPrime[MAX];
vector<int> Era(int n = MAX) {
    vector<int> res;
    IsPrime[0] = false; IsPrime[1] = false;
    for (int i = 2; i < n; ++i) IsPrime[i] = true;
    for (int i = 2; i < n; ++i) {
        if (IsPrime[i]) {
            res.push_back(i);
            for (int j = i*2; j < n; j += i) IsPrime[j] = false;
        }
    }
    return res;
}

int N;
vector<long long> a;
vector<int> prs;

long long calc(long long n) {
    long long res = 1;
    for (auto p : prs) {
        int num = 0;
        while (n % p == 0) ++num, n /= p;
        if (num == 0) continue;
        else if (num == 1) res *= p * p;
        else res *= p;
    }
    long long left = 1, right = 1000000;
    while (right - left > 1) {
        long long mid = (left + right) / 2;
        if (mid * mid > n) right = mid;
        else left = mid;
    }
    if (left * left == n) {
        res *= left;
        return res;
    }
    else {
        if (n > 1000000) return -1;
        res *= n * n;
        if (res > 10000000000LL) return -1;
        return res;
    }
}

long long solve() {
    int res = 0;
    prs = Era();
    map<long long, int> ma;
    for (int i = 0; i < N; ++i) {
        long long tmp = 1;
        for (auto p : prs) {
            int num = 0;
            while (a[i] % p == 0) ++num, a[i] /= p;
            num %= 3;
            for (int j = 0; j < num; ++j) tmp *= p;
        }
        tmp *= a[i];
        
        if (tmp == 1) res = 1;
        else ma[tmp]++;
    }
    for (auto it = ma.begin(); it != ma.end(); ++it) {
        long long val = it->first;
        long long aite = calc(val);
        int tmp = max(it->second, ma[aite]);
        it->second = ma[aite] = 0;
        res += tmp;
    }
    return res;
}

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