题目
题目背景
非洲猪瘟爆发,村里仅剩的小香猪肉价疯涨不下!以前经常买肉的 O n e I n D a r k \sf OneInDark OneInDark 掏空了腰包,还是只能偶尔揩揩油……
题目描述
这只小香猪非常狡猾,它要利用自己的身姿,牟取最多的利润!
现在有 2 n 2n 2n 个人,都想来买肉。小香猪决定选出其中 n n n 个人,但是他们付的钱会是平常肉价的 gcd λ i \gcd\lambda_i gcdλi 倍,其中 λ i \lambda_i λi 是第 i i i 个人的男性化程度(因为小香猪喜欢当妹妹)。
那么小香猪可以做到的最大倍率是多少呢?
数据范围与提示
n ⩽ 1 0 5 n\leqslant 10^5 n⩽105 但 1 ⩽ λ i ⩽ 1 0 9 1\leqslant\lambda_i\leqslant 10^9 1⩽λi⩽109 。
思路
最关键的就是 2 n 2n 2n 中选出 n n n 个。这有什么用呢?
最初以为是鸽笼原理:除去 ( n + 1 ) (n+1) (n+1) 种特殊情况,总是存在两个相邻的人同时被选中。然而这并没有什么用处……
想了很久,灵机一动想起了 随机化。随机一个人 δ \delta δ,它在最优方案中存在的概率是 1 2 \frac{1}{2} 21,多随机几次,就可以保证找到一个最优方案中的人。于是只用考虑 δ \delta δ 的因数是否能作为答案。
最暴力的方法是,对于 d d d,找找有多少个数是 d d d 的倍数。而我们检测的所有 d d d 都是 δ \delta δ 的因数,所以先将其余的数与 δ \delta δ 取 gcd \gcd gcd 不会影响判断。那么取完 gcd \gcd gcd 之后,得到的结果都是 δ \delta δ 的因数,问题转化为,给出 δ \delta δ 若干因数,求每个因数的所有倍数的个数之和。
如果将所有数质因数分解,将每个质因子视作一个维度,其指数是坐标;那么相当于求一个高维偏序,做后缀和即可。由于只有 log A \log A logA 维,共 A \sqrt{A} A 个因子,时间复杂度 O ( A log A ) \mathcal O(\sqrt{A}\log A) O(AlogA) 。
当然,我尝试了两种写法:将所有因数丢进 m a p \tt map map,然后 2 n 2n 2n 个 gcd \gcd gcd 的结果就直接映射到一个高维坐标(用一个整数状压表示);或者将 2 n 2n 2n 个 gcd \gcd gcd 的结果拿去做质因数分解。按理说,前者应该是 O ( n log A + A log A ) \mathcal O(n\log A+\sqrt{A}\log A) O(nlogA+AlogA) 的,而后者是 O ( n log A + A ) \mathcal O(n\log A+\sqrt{A}) O(nlogA+A) 的(如果不考虑求后缀和的部分)。可实地测试中,前者更快!
原因大概是因为,因数个数其实比 A \sqrt A A 小。我之前看这里面说大概是 O ( A ) \mathcal O(\sqrt A) O(A) 的,现在再看,发现它还带一个 1 10 1\over 10 101 的系数……果然大 O \mathcal O O 就是大 O \mathcal O O,从来都不像 Θ \Theta Θ 那样可信……
代码
#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cctype>
#include <random>
#include <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 llong getGcd(llong a,llong b){
for(; true; b%=a){
if(b == 0) return a;
if((a %= b) == 0) return b;
}
return 0; // what the heck
}
const int MAXN = 200005;
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]<=n/i; ++j)
isPrime[i*primes[j]] = false;
}
}
const int LOGA = 50;
int zs[LOGA], coe[LOGA], tot;
map<llong,int> haxi; ///< from divisor to index
void get_factor(llong x){
tot = 0; coe[0] = 1, zs[0] = 0;
static vector< pair<llong,int> > _tmp;
haxi.clear(); _tmp.clear();
haxi[1] = 0; // the only survivor
for(int i=1; i<=primes_size&&primes[i]<=x/primes[i]; ++i){
if(x%primes[i]) continue; // no such thing
for(++tot,zs[tot]=0; !(x%primes[i]); )
x /= primes[i], ++ zs[tot];
coe[tot] = coe[tot-1]*(zs[tot-1]+1);
for(auto it : haxi){
llong v = it.first; int _id = it.second;
for(int j=1; j<=zs[tot]; ++j){
v *= primes[i], _id += coe[tot];
_tmp.push_back(make_pair(v,_id));
}
}
for(; !_tmp.empty(); _tmp.pop_back())
haxi[_tmp.back().first] = _tmp.back().second;
}
if(x == 1) return ;
zs[++ tot] = 1, coe[tot] = coe[tot-1]*(zs[tot-1]+1);
for(auto it : haxi) // consider this big divisor
_tmp.push_back(make_pair(it.first*x,it.second+coe[tot]));
for(; !_tmp.empty(); _tmp.pop_back())
haxi[_tmp.back().first] = _tmp.back().second;
}
int cnt[SQRTA<<1];
int main(){
sieve(); // maybe optimization backwards ...
mt19937 rnd; // rander
rnd.seed((unsigned)1983);
int n = int(readint());
uniform_int_distribution<int> _sy(1,n<<1);
rep(i,1,n<<1) a[i] = readint();
llong ans = 1; // real answer
for(int cs=20; cs; --cs){
const int pivot = _sy(rnd);
get_factor(a[pivot]);
const int all = coe[tot+1] = coe[tot]*(zs[tot]+1);
memset(cnt,0,all<<2); // clear
rep(i,1,n<<1) ++ cnt[haxi[getGcd(a[i],a[pivot])]];
rep(i,1,tot) drep(j,all-1,0)
if((j%coe[i+1])/coe[i] != zs[i])
cnt[j] += cnt[j+coe[i]];
for(auto it=haxi.rbegin(); it!=haxi.rend(); ++it)
if(cnt[it->second] >= n){
ans = max(ans,it->first); break;
}
}
printf("%lld\n",ans);
return 0;
}