2017-2018 Petrozavodsk Winter Training Camp, Saratov SU Contest F. GCD

問題文

http://codeforces.com/gym/101741/problem/F

問題概要

長さ n の数列のうち、k個を削除して残ったものの GCD の最大値を求めよ.

• n <= 10^5
• k <= n / 2
• a_i <= 10^18

解法

まず, 残すもののうちの１つを乱択で決め打ちする.
k が n / 2 以下であるという制約からこれは確率 1/2 なので10回くらい試せばある程度安心.

こうして１つを決めうちする (その値を X とする) と, 他の要素は X との GCD であると見なせるので, Xの約数のみを考えればよくなる.

とりあえず X を素因数分解する必要があるが, 愚直にやると間に合わないので ミラーラビン法 + ポラード・ロー法 を使うと O(X^(1/4) logX) くらいで求まる.

あとは, 倍数のゼータ変換 (?) を行えばよくて, これは各素因数ごとに dfs して大きい方から足し合わせていけば O( (Xの約数の個数) * (X の素因数の種類数) ) くらいになる.
大きそうに見えるが, 10^18 程度の値なら約数の個数は案外小さいので, なんとか間に合う.

ソースコード

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

ll mod_prod(ll a, ll b, ll md) {
    ll res = (a * b - (ll)((long double)a / md * b) * md) % md;
    return res < 0 ? res + md : res;
}

ll mod_pow(ll x, ll k, ll m) {
    ll res = 1;
    while (k) {
        if (k & 1) res = mod_prod(res, x, m);
        k /= 2;
        x = mod_prod(x, x, m);
    }
    return res;
}

bool miller_rabin_test(ll n) {
    const static vector<ll> sprp = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 };
    if (n < 2) return false;
    ll d = n - 1, s = 0;
    while ((d & 1) == 0) d >>= 1, ++s;
    for (ll a : sprp) if (a <= n) {
        if (a == n) return true;
        if (mod_pow(a, d, n) != 1) {
            bool ok = true;
            for (ll r = 0; r < s; r++) {
                if (mod_pow(a, d << r, n) == n - 1) {
                    ok = false;
                    break;
                }
            }
            if (ok) return false;
        }
    }
    return true;
}

ll pollard_rho(ll n) {
    static mt19937 mt(random_device{}());
    if (n == 4) return 2;
    ll c = mt() % n;
    ll x = mt() % n, y = x, d = 1;
    while (d == 1) {
        x = (mod_prod(x, x, n) + c) % n;
        y = (mod_prod(y, y, n) + c) % n;
        y = (mod_prod(y, y, n) + c) % n;
        d = __gcd(abs(x - y), n);
    }
    return d;
}

template <typename T>
vector<pair<T, int>> factorize(T n) {
    if (n <= 1) return{};
    vector<T> st = { n }, ps;
    while (!st.empty()) {
        T t = st.back(); st.pop_back();
        if (t <= 1) continue;
        if (miller_rabin_test(t)) {
            ps.push_back(t);
            continue;
        }
        T x;
        while (x = pollard_rho(t), x == t);
        st.push_back(x);
        st.push_back(t / x);
    }
    sort(ps.begin(), ps.end());
    vector<pair<T, int>> res;
    for (auto p : ps) {
        if (res.empty() || res.back().first != p) res.emplace_back(p, 1);
        else res.back().second++;
    }
    return res;
}

mt19937 mt(random_device{}());

int main()
{
    int n, k;
    cin >> n >> k;
    vector<ll> a(n);
    for (int i = 0; i < n; i++) {
        cin >> a[i];
    }
    ll res = 0;
    vector<int> used(n);
    for (int cc = 0; cc < min(n, 10); cc++) {
        int id;
        while (used[id = mt() % n]);
        used[id] = 1;
        unordered_map<ll, int> cnt;
        for (int i = 0; i < n; i++) {
            cnt[__gcd(a[i], a[id])] += 1;
        }
        auto fs = factorize(a[id]);
        int sz = fs.size();
        for (auto pp : fs) {
            ll p = pp.first;
            function<void(int, ll)> calc = [&](int it, ll val) {
                if (it == sz) {
                    for (int i = 0; i < pp.second; i++) {
                        val *= p;
                    }
                    for (int i = pp.second - 1; i >= 0; i--) {
                        cnt[val / p] += cnt[val];
                        if (i > 0) val /= p;
                    }
                    return;
                }
                if (fs[it].first == p) {
                    calc(it + 1, val);
                    return;
                }
                for (int i = 0; i <= fs[it].second; i++) {
                    calc(it + 1, val);
                    if (i + 1 <= fs[it].second) val *= fs[it].first;
                }
            };
            calc(0, 1);
        }
        for (auto p : cnt) if (p.second >= n - k) {
            res = max(res, p.first);
        }
    }
    cout << res << endl;
    return 0;
}

感想

久々にためになる問題だと思った.
10^18 くらいの値は定数個くらいなら素因数分解できるのとか, 約数の個数がある程度小さいのとかは典型っぽい気がする.