CF364D Ghd 题解

题目链接：CF364D Ghd

题意：

定义一组数 \(a_1, a_2, \cdots, a_n\) 的半最大公约数 (Greatest half-common divisor, GHD) 为最大的整数 \(d\)，满足 \(a_1, a_2, \cdots, a_n\) 中至少有 \(\left \lceil \frac{n}{2} \right \rceil\) 个 \(d\) 的倍数。

给定 \(n\) 个正整数 \(a_1, a_2, \cdots, a_n\)，请你求出它们的半最大公约数。

输入格式：

第一行包含一个正整数 \(n~(1\le n \leq 10^6)\)，表示数的个数。

第二行包含 \(n\) 个正整数 \(a_1, a_2, \cdots, a_n(1 \leq a_i \leq 10^{12})\)，表示这 \(n\) 个数的值。

输出格式：

输出一行一个整数，表示这些数的半最大公约数。

前置知识：\(d(10^{12}) \le 6720\) 。

注意到最终的 GHD 的集合 \(S\) 中，一定存在至少一个 \(M\) ，使得 \(\forall x \in S\) 有 \(x \mid M\) 。

因此最简单的思路就是固定一个数 \(M = a_k\) ，然后枚举 \(g \mid M\) ，看看有多少个 \(a_i\) 满足 \(g \mid a_i\) 。

如果个数大于等于 \(\left\lceil\frac{n}{2}\right\rceil\) ，那么我们就 \(\mathtt{res}\uparrow g\) 一下。这样做的复杂度是 \(\mathcal{O}(d(M)\times n^2)\) 的。

考虑优化这个做法。注意到当 \(g \mid M\) 时，\(g \mid a_i\) 当且仅当 \(g \mid \gcd(M,a_i)\) 。

因此我们可以令 \(a_i \gets \gcd(M,a_i)\) ，这样所有的 \(a_i\) 都是 \(M\) 的约数了。这个时候就相对好办些了。

记 \(d_i\) 表示 \(M\) 的第 \(i\) 个约数，并记 \(f_i\) 表示有多少个 \(a_i=d_i\) （ \(a_i\) 是修改后的，且 \((0\le i <d(M)\) ）

容易发现，此时的 \(\sum [g\mid a_i]\) 就可以用下面的式子快速得到（不妨设 \(g=d_i\) ） \[ \sum_{i \le j < d(M)}f_j \times [d_i \mid d_j] \] 记 \(m = d(M)\) ，则此时的复杂度为 \(\mathcal{O}\left(nm^2 + n^2 \log m\right)\) ，复杂度瓶颈仍在于枚举 \(M\) 。

不过到现在为止，我们还没有用到 \(|S| \geq\left\lceil\frac{n}{2}\right\rceil\) 这个特殊性质。

可以发现在这个条件下，任意一个元素 \(a_i \in S\) 的概率不小于 \(\frac{1}{2}\) 。

因此我们可以随机选择一个 \(M\) ，设随机 \(k\) 次，则算法的正确率为（考虑每个 \(M\) 都不在最终 GHD 时的概率） \[ 1 - \frac{1}{2^k} \] 时间复杂度 \(\mathcal{O}\left(k(m^2 + n \log m)\right)\) ，\(k=12\) 的时候差不多。

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(1e6+15))
#define D ((int)(7015))

const int k = 12;
int n,res = 1,f[D],a[N],d[D];
int gcd(int a,int b) { return b == 0 ? a : gcd(b, a % b); }
int solve(int x)
{
    int cnt = 0, sq = sqrt(x);
    for(int i=1; i * i < x; i++)
        if(x % i == 0) { d[cnt++] = i; d[cnt++] = x/i; }
    if(sq * sq == x) d[cnt++] = sq;
    sort(d, d+cnt); for(int i=0; i<cnt; i++) f[i] = 0;
    for(int i=0; i<n; i++) ++f[lower_bound(d, d+cnt, gcd(a[i],x)) - d];
    for(int i=cnt-1,tot = 0; i; i--, tot = 0)
    {
        for(int j=i; j<cnt; j++) if(d[j] % d[i] == 0)  tot += f[j];
        if(tot >= (n + 1) / 2) return d[i];
    }
    return 1;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    char *_ = new char;
    mt19937 rd(time(0) + (unsigned int)_); delete _;
    cin >> n; for(int i=0; i<n; i++) cin >> a[i];
    for(int i=0; i<k; i++) up(res, solve(a[rd() % n]));
    cout << res << '\n';
    return 0;
}