概要

この問題は、与えられた \(N\) 個の宝石から \(2\) 個以上を選び、それらの品質値の最大公約数（GCD）を最大化し、さらにその最大GCDを達成するために必要な最小の宝石の個数 \(K\) を求め、使わなかった宝石の個数 \(N - K\) を出力する問題です。

品質値の最大値を \(M\) とすると、最大公約数の最大値 \(G\) の決定、および包除原理を用いた \(K\) の探索を、モジュロハッシュを組み合わせることで効率的に行うことができます。

---

考察

1. 最大GCD \(G\) の決定

選んだ宝石たちのGCDが \(g\) になるためには、選んだすべての宝石の品質値が \(g\) の倍数 である必要があります。 したがって、「\(A\) の中に \(g\) の倍数が \(2\) 個以上存在する」ような最大の \(g\) が、達成可能な最大の美しさ \(G\) になります。

これは、品質値の最大値 \(M = \max(A_i)\) から降順に \(g = M, M-1, \ldots, 1\) と走査し、各 \(g\) について「\(A\) に含まれる \(g\) の倍数の個数」をカウントすることで求められます。最初に対象の個数が \(2\) 以上になった \(g\) が求める \(G\) です。

2. 問題の単純化（GCDが1となる最小個数 \(K\) の探索）

最大GCD \(G\) が決まれば、品質値が \(G\) の倍数である宝石のみを考慮すればよいことになります。 これらの宝石の品質値を \(G\) で割った値を \(C_1, C_2, \ldots, C_m\) とします。

元の問題は、以下のように言い換えることができます。

「数列 \(C\) から \(k\) 個の要素を選んで、そのGCDが \(1\) になるような最小の \(k\)（これを \(K\) とする）を求める」

なぜなら、割った後のGCDが \(1\) になるということは、元の品質値のGCDがちょうど \(G\) になるということだからです。

具体例

\(A = [12, 18, 30]\) の場合： 1. \(g=6\) の倍数は \([12, 18, 30]\) の \(3\) 個（\(2\) 個以上）存在するため、最大GCDは \(G = 6\) となります。 2. これらを \(G=6\) で割ると、 \(C = [2, 3, 5]\) になります。 3. \(C\) の中から \(k\) 個選んでGCDを \(1\) にする最小の個数を考えます。 - \(k=2\) のとき、例えば \(\{2, 3\}\) を選ぶと \(\gcd(2, 3) = 1\) となり、GCDを \(1\) にできます。 - したがって、最小個数 \(K = 2\) となり、答えは \(N - K = 3 - 2 = 1\) です。

3. GCDが1になる最小個数 \(K\) の高速な判定

\(k = 2, 3, \ldots\) と増やしていき、「\(C\) から \(k\) 個選んでGCDが \(1\) になる組み合わせが \(1\) 通り以上存在するか」を判定します。

「GCDがちょうど \(1\) になる選び方の数」を直接数えるのは難しいですが、包除原理（メビウスの反転公式） を用いることで高速に計算できます。

\(C\) の中で \(x\) の倍数であるものの個数を \(cnt[x]\) とします。 \(C\) から \(k\) 個選んで、そのGCDが \(x\) の倍数であるような組み合わせの数は \(\binom{cnt[x]}{k}\) 通りです。 これにメビウス関数 \(\mu(x)\) を掛け合わせて足し合わせることで、選んだ \(k\) 個のGCDがちょうど \(1\) になる組み合わせの総数 \(S_k\) は次のように求まります。

\[S_k = \sum_{x=1}^{V} \mu(x) \binom{cnt[x]}{k}\]

（ここで \(V = \max(C)\)）

\(S_k > 0\) となる最小の \(k\) が求める \(K\) です。 ただし、 \(S_k\) は非常に大きな値になるため、適切な大きな素数 \(P_1, P_2\) を用いたモジュロ演算（ダブルハッシュ）を行い、 \(S_k \not\equiv 0 \pmod{P_1}\) または \(S_k \not\equiv 0 \pmod{P_2}\) であるかを確認することで、ハッシュ衝突を防ぎながら正しく判定できます。

---

アルゴリズム

1. 頻度配列の作成: 各品質値 \(A_i\) の出現回数を記録する配列 count_A を作成します。
2. 最大GCD \(G\) の探索: \(g\) を \(M\) から \(1\) まで降順にループし、調和級数の要領で \(g\) の倍数の個数を足し合わせます。個数が \(2\) 以上になった最初の \(g\) を \(G\) とします。
3. 数列 \(C\) の構築: \(A_i\) のうち \(G\) の倍数であるものを \(G\) で割り、数列 \(C\) を作ります。 \(V = \max(C)\) とします。
4. 倍数カウンタ cnt の作成: \(C\) の要素について、各 \(x \ (1 \le x \le V)\) の倍数が何個あるかを調和級数を用いて \(O(V \log V)\) で求め、 cnt[x] に格納します。
5. メビウス関数の計算: 線形篩（Linear Sieve）を用いて、 \(1\) から \(V\) までのメビウス関数 \(\mu(x)\) を \(O(V)\) で事前計算します。
6. 最小個数 \(K\) の探索: \(k = 2, 3, \ldots\) について、 \(S_k \pmod{P_1}\) および \(S_k \pmod{P_2}\) を計算します。 どちらか一方でも \(0\) でなければ、GCDが \(1\) になる \(k\) 個の選び方が存在することになるため、その \(k\) を \(K\) として探索を終了します。

【ポイント】 \(K\) の探索範囲について、任意の \(10^6\) 以下の整数が持つ異なる素因数の個数は高々 \(7\) 個（\(2 \times 3 \times 5 \times 7 \times 11 \times 13 \times 17 > 10^6\) より）です。したがって、互いに素なグループを作るために必要な要素数 \(K\) は非常に小さな値（高々 \(7\) 程度）に収まります。そのため、 \(k\) のループは極めて少ない回数で終了します。

---

計算量

\(M = \max(A_i) \le 10^6\) とします。

• 時間計算量: \(O(M \log M)\)

  • \(G\) の決定に \(O(M \log M)\)（調和級数の和 \(\sum \frac{M}{g} \approx M \log M\)）
  • cnt 配列の計算に \(O(V \log V)\)
  • メビウス関数の計算に \(O(V)\)
  • \(K\) の探索ループは、各ステップで \(O(V)\) の計算を行います。ループ回数 \(K\) は定数（最大でも \(7\) 程度）とみなせるため、この部分は \(O(V)\) となり、全体として \(O(M \log M)\) で十分に高速に動作します。

• 空間計算量: \(O(M)\)

  • 値の頻度やメビウス関数を保持する配列のサイズは最大で \(M\) であるため、メモリ制限に対しても非常に余裕があります。

---

実装のポイント

• ダブルハッシュによる衝突回避: 組み合わせ数 \(S_k\) が非常に大きくなるため、 \(1\) つの素数によるモジュロ演算だけでは、たまたま \(S_k \equiv 0 \pmod P\) となってしまう「偽陰性（本来は存在するのに存在しないと判定される）」が発生する可能性があります。 \(10^9+7\) と \(10^9+9\) のような \(2\) つの異なる素数でダブルハッシュを取ることで、衝突の確率を実質的にゼロにできます。

• 高速な二項係数 \(\binom{n}{k} \pmod P\) の計算: \(k\) が小さいため、分母の階乗の逆元（inv_k_fact）をあらかじめ求めておくことで、分子を \(O(k)\) で掛け合わせるだけで高速に \(\binom{n}{k}\) を計算できます。

  ソースコード

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

long long power(long long base, long long exp, long long mod) {
    long long res = 1;
    base %= mod;
    while (exp > 0) {
        if (exp % 2 == 1) res = (res * base) % mod;
        base = (base * base) % mod;
        exp /= 2;
    }
    return res;
}

long long modInverse(long long n, long long mod) {
    return power(n, mod - 2, mod);
}

long long nCr(long long n, int k, long long inv_k_fact, long long mod) {
    if (n < k) return 0;
    long long num = 1;
    for (int i = 0; i < k; ++i) {
        num = (num * ((n - i) % mod)) % mod;
    }
    return (num * inv_k_fact) % mod;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int N;
    if (!(cin >> N)) return 0;

    vector<int> A(N);
    int max_A = 0;
    for (int i = 0; i < N; ++i) {
        cin >> A[i];
        if (A[i] > max_A) max_A = A[i];
    }

    vector<int> count_A(max_A + 1, 0);
    for (int x : A) {
        count_A[x]++;
    }

    int G = 1;
    for (int g = max_A; g >= 1; --g) {
        int multiples = 0;
        for (int y = g; y <= max_A; y += g) {
            multiples += count_A[y];
        }
        if (multiples >= 2) {
            G = g;
            break;
        }
    }

    vector<int> C;
    for (int x : A) {
        if (x % G == 0) {
            C.push_back(x / G);
        }
    }

    int V = 0;
    for (int x : C) {
        if (x > V) V = x;
    }

    vector<int> count_C(V + 1, 0);
    for (int x : C) {
        count_C[x]++;
    }

    vector<int> cnt(V + 1, 0);
    for (int x = 1; x <= V; ++x) {
        for (int y = x; y <= V; y += x) {
            cnt[x] += count_C[y];
        }
    }

    vector<int> mu(V + 1, 0);
    vector<int> primes;
    vector<bool> is_prime(V + 1, true);
    mu[1] = 1;
    for (int i = 2; i <= V; ++i) {
        if (is_prime[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            if (i * p > V) break;
            is_prime[i * p] = false;
            if (i % p == 0) {
                mu[i * p] = 0;
                break;
            } else {
                mu[i * p] = -mu[i];
            }
        }
    }

    long long P1 = 1000000007;
    long long P2 = 1000000009;

    int K = -1;
    for (int k = 2; k <= (int)C.size(); ++k) {
        long long fact1 = 1, fact2 = 1;
        for (int i = 1; i <= k; ++i) {
            fact1 = (fact1 * i) % P1;
            fact2 = (fact2 * i) % P2;
        }
        long long inv_k_fact1 = modInverse(fact1, P1);
        long long inv_k_fact2 = modInverse(fact2, P2);

        long long sum1 = 0;
        long long sum2 = 0;
        for (int x = 1; x <= V; ++x) {
            if (mu[x] == 0) continue;
            long long term1 = nCr(cnt[x], k, inv_k_fact1, P1);
            long long term2 = nCr(cnt[x], k, inv_k_fact2, P2);
            if (mu[x] == 1) {
                sum1 = (sum1 + term1) % P1;
                sum2 = (sum2 + term2) % P2;
            } else {
                sum1 = (sum1 - term1 + P1) % P1;
                sum2 = (sum2 - term2 + P2) % P2;
            }
        }
        if (sum1 != 0 || sum2 != 0) {
            K = k;
            break;
        }
    }

    cout << N - K << "\n";

    return 0;
}

---

この解説は gemini-3.5-flash-thinking によって生成されました。