[CF364D]Ghd

题目

传送门 to CF：下为等价题面。

题目背景
非洲猪瘟爆发，村里仅剩的小香猪肉价疯涨不下！以前经常买肉的 $\mathsf{O}\mathsf{n}\mathsf{e}\mathsf{I}\mathsf{n}\mathsf{D}\mathsf{a}\mathsf{r}\mathsf{k}$ 掏空了腰包，还是只能偶尔揩揩油……

题目描述
这只小香猪非常狡猾，它要利用自己的身姿，牟取最多的利润！

现在有 $n$ 个人，都想来买肉。小香猪决定选出其中 $\lceil \frac{n}{2}\rceil$ 个人，但是他们付的钱会是平常肉价的 $\gcd {\lambda }_i$ 倍，其中 ${\lambda }_i$ 是第 $i$ 个人的男性化程度（因为小香猪喜欢当妹妹）。

那么小香猪可以做到的最大倍率是多少呢？

数据范围与提示
$n⩽10^6$ 但 $1⩽{\lambda }_i⩽10^{12}$ 。

思路

最关键的就是 $n$ 中选出 $\frac{n}{2}$ 个。这有什么用呢？

最初以为是鸽笼原理：除去一些特殊情况，总是存在两个相邻的人同时被选中。然而这并没有什么用……

想了很久，灵机一动想起了 随机化。随机一个人 $\delta$，它在最优方案中存在的概率是 $\frac{1}{2}$，多随机几次，就可以保证找到一个最优方案中的人。于是只用考虑 $\delta$ 的因数是否能作为答案。

最暴力的方法是，对于 $d$，找找有多少个数是 $d$ 的倍数。而我们检测的所有 $d$ 都是 $\delta$ 的因数，所以先将其余的数与 $\delta$ 取 $\gcd$ 不会影响判断。那么取完 $\gcd$ 之后，得到的结果都是 $\delta$ 的因数，问题转化为，给出 $\delta$ 若干因数，求每个因数的所有倍数的个数之和。

如果将所有数质因数分解，将每个质因子视作一个维度，其指数是坐标；那么相当于求一个高维偏序，做后缀和即可。由于只有 $\log A$ 维，共 $\sqrt{A}$ 个因子，时间复杂度 $\mathcal{O}(\sqrt{A}\log A)$ 。

当然，因数个数其实比 $\sqrt{A}$ 小。我之前看这里面说大概是 $\mathcal{O}(\sqrt{A})$ 的，现在再看，发现它还带一个 $\frac{1}{10}$ 的系数……

最离谱的是，$\text{CF}$ 上只能随机化 $10$ 次，不然就会 $\mathtt{T}\mathtt{L}\mathtt{E}$ 了！而且，许多 $\mathcal{O}(d^2)$ 检查的也能通过！其中 $d$ 为因数数量。我很难理解！

代码

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cctype>
#include <random>
#include <unordered_map>
using namespace std;
#define rep(i,a,b) for(int i=(a); i<=(b); ++i)
#define drep(i,a,b) for(int i=(a); i>=(b); --i)
typedef long long llong;
inline llong readint(){
	llong a = 0; int c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar())
		if(c == '-') f = -f;
	for(; isdigit(c); c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}
inline void writeint(const llong &x){
	if(x > 9) writeint(x/10);
	putchar(int((x%10)^48));
}
inline llong getGcd(llong a,llong b){
	for(; b; a%=b,swap(a,b));
	return a; // when b == 0
}

const int MAXN = 1000005;
llong a[MAXN];

const int SQRTA = 1000000;
int primes[SQRTA+5], primes_size;
bool isPrime[SQRTA+5];
void sieve(int n = SQRTA){
	memset(isPrime+2,true,n-1);
	for(int i=2; i<=n; ++i){
		if(isPrime[i]) primes[++ primes_size] = i;
		for(int j=1; j<=primes_size&&primes[j]*i<=n; ++j){
			isPrime[i*primes[j]] = false;
			if(i%primes[j] == 0) break;
		}
	}
}

const int LOGA = 50;
llong fac[SQRTA], pri_fac[LOGA]; int tot, cnt;
void get_factor(llong x){
	static llong _tmp[SQRTA];
	int _len = 1; _tmp[1] = x;
	llong y = x; tot = cnt = 0;
	for(int i=2; llong(i)*i<=x; ++i) if(!(x%i)){
		fac[++ cnt] = i; // divisor
		if(llong(i)*i != x) _tmp[++ _len] = x/i;
		if(!isPrime[i] || y%i != 0) continue;
		pri_fac[++ tot] = i, y /= i;
		while(!(y%i)) y /= i; // get primes
	}
	if(y != 1) pri_fac[++ tot] = y;
	rep(i,1,_len) fac[cnt+i] = _tmp[_len+1-i];
	std::inplace_merge(fac+1,fac+_len+1,
		fac+_len+cnt+1); cnt += _len;
}

std::unordered_map<llong,int> mp;
llong ans; int hit[SQRTA];
int main(){
	sieve(); std::mt19937 rnd;
	rnd.seed((unsigned)14787841);
	int n = int(readint());
	uniform_int_distribution<int> _sy(1,n);
	rep(i,1,n) a[i] = readint();
	ans = 1; // forever answer
	for(int cs=10; cs; --cs){
		const int pivot = _sy(rnd);
		get_factor(a[pivot]); // core
		mp.clear(); mp.reserve(cnt+1);
		rep(i,1,cnt) mp[fac[i]] = i;
		mp[1] = 0; memset(hit+1,0,cnt<<2);
		rep(i,1,n) ++ hit[mp[getGcd(a[i],a[pivot])]];
		llong *end_j = pri_fac+tot+1;
		for(llong *j=pri_fac+1; j!=end_j; ++j)
			for(int i=cnt; fac[i]>ans; --i)
				if(!((a[pivot]/fac[i])%(*j)))
					hit[i] += hit[mp[(*j)*fac[i]]];
		for(int i=cnt; fac[i]>ans; --i)
			if((hit[i]<<1) >= n){
				ans = fac[i]; break;
			}
	}
	writeint(ans), putchar('\n');
	return 0;
}