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