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

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

Codeforces #613 (Div. 2) F. Classical? (R2800)

勉強になった...けど、これ知らずにできるもんなの!?

問題へのリンク

あと、LCM の最小値バージョンもある!

drken1215.hatenablog.com

問題概要

 N 個の正の整数  a_{1}, a_{2}, \dots, a_{N} が与えられる。これらから 2 個選んで LCM をとってできる  \frac{N(N-1)}{2} 個の整数の最大値を求めよ。

制約

  •  2 \le N \le 10^{5}
  •  1 \le a_{i} \le 10^{5}

考えたこと

 a_{i} \le 10^{5} という制約がいかにもそれっぽい。LCM ってのはなかなかに考えづらいから、ここは GCD を固定して考えることにする。つまり

  • GCD の値が  g となるような 2 数についての LCM の最大値

を各  g について求めていくことにする。 g \le 10^{5} なのでいい感じ。さて、考察対象を  g の倍数に絞ったとき、それらの中から「互いに素」な 2 つの整数の積として考えられる最大値を求めたい。

まずこの手の問題でよくあることとして、数列  aバケットで管理しておくことにする。つまり

  • a'[  v ] := 数列  a の中に  v という値があるかどうか

という配列として管理する。こうすることで、 g を抜き取る作業を  O(A/g) でできるのだ ( A を登場する値の最大値とする)。 g = 1, 2, \dots, A で合計したとき、 O(A\log{A}) となる。

g の倍数に関する問題

問題は  m = O(A/g) 個の値  b_{0}, b_{1}, \dots, b_{m-1} から互いに素な 2 整数の積の最大値を考える問題に帰着された。

注目したいこととして、

  •  a, b が互いに素であったとき、 b の相方として  a より小さいものは考えなくてよい

ということ。当たり前のことなのだが、たったこれだけでかなり枝刈りができる。以下のようにする。


  • 大きい順に値を見ていってスタックに入れる
  • 今見ている値が  x であるとき、スタックの中に  x と互いに素な値があるならば、スタックの中で最小の値を破棄してよい

その理由は、もしスタックの中に値  c が存在してそれが  x と互いに素であるならば、

  •  c から見て  x より小さい値を考慮する必要がないため、 c のことは今後一切考えなくてもよいので  c を捨てることができる
  • スタックにある値のうち  c より小さい値については、今後  x より小さい値としか結びつく可能性がないため、同様に捨てることができる

という感じ。結局、どの要素もスタックに 1 回だけ入れられて、1 回だけ削除されるような感じになる。すごくこどふぉっぽい計算量評価の仕組みだ!!!!!!!

残る問題は「スタックの中に  x と互いに素な整数があるかどうか」をどうやって高速に判定するかだ。そういうことができるデータ構造を設計しよう。

やりたいこと

やりたいことは以下の操作のできるデータ構造を設計すること。登場する値の最大値は  A 程度であるとする。

  1.  x を挿入する
  2.  x を削除する
  3. データ構造の中に値  x と互いに素なものが何個あるかを求める (存在判定だけでよいけど、個数も求められる)

これを実現するために、以下のデータを管理する。

  • counter[  v ] := データ構造の中に  v の倍数が何個あるか

 x の挿入は、 x のすべての約数  d に対して counter[  d ] をインクリメントすればよい。 x の削除は、 x のすべての約数  d に対して counter[  d ] をデクリメントすればよい。最後に  x と互いに素なものの個数は、

 \sum_{d | x} \mu(d) \times {\rm counter} \lbrack d \rbrack

で求められる。 \mu(d)メビウス関数。そしてこれらのクエリを処理するためには

を前処理しておく。これらはエラトステネスの篩によってできる。

計算量

整数  1, \dots, A の約数の個数の最大値を  d(A) とする。結局計算量は、「約数の個数分」だけの時間のかかる処理を  O(A\log{A}) 回やる感じになるので、 O(A d(A) \log{A}) みたいな感じになると思う。

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

// 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;
    }
};



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

const int MAX = 110000;

int main() {
    Eratos er(MAX);
    vector<vector<int>> divs(MAX);
    for (int n = 1; n < MAX; ++n) divs[n] = er.divisors(n);

    stack<int> st;
    vector<int> counter(MAX, 0);
    auto push = [&](int x) {
        st.push(x);
        for (auto d : divs[x]) ++counter[d];
    };
    auto pop = [&]() {
        int x = st.top();
        st.pop();
        for (auto d : divs[x]) --counter[d];
    };
    auto count = [&](int x) {
        int res = 0;
        for (auto d : divs[x]) res += er.mebius[d] * counter[d];
        return res;
    };

    int N;
    cin >> N;
    vector<long long> a(N);
    vector<bool> isin(MAX+1, false);
    long long res = 0;
    for (int i = 0; i < N; ++i) {
        cin >> a[i];
        res = max(res, a[i]);
        isin[a[i]] = true;
    }

    for (int g = 1; g < MAX; ++g) {
        for (int v = (MAX/g)*g; v >= g; v -= g) {
            if (!isin[v]) continue;
            while (count(v/g)) {
                long long t = st.top() * g;
                res = max(res, t / GCD(v, t) * v);
                pop();
            }
            push(v/g);
        }
        while (!st.empty()) pop();
    }
    cout << res << endl;
}