给定 n(2≤n≤5e5)n(2\le n\le 5e5)n(2≤n≤5e5) 个数,2≤ai≤1e72\le a_i\le 1e72≤ai≤1e7。
一次操作中,从数组 aaa 中选定一个 xxx,再从剩下的数中选出一个 gcd>1\gcd>1gcd>1 的非空集合 YYY,使得所有选中的数的 gcd=1\gcd=1gcd=1,即 gcd(x,Y)=1\gcd(x,Y)=1gcd(x,Y)=1。
问有多少种这样的方案,答案 mod 1e9+7mod\ 1e9+7mod 1e9+7。
我们尝试用公式表达出我们的所求:
∑Y⊆U∑x⊈Y[gcd(Y)≠1][gcd(x,Y)=1]\sum_{Y \subseteq U}\sum_{x\not\subseteq Y}[\gcd(Y)\neq1][\gcd(x,Y)=1]Y⊆U∑x⊆Y∑[gcd(Y)=1][gcd(x,Y)=1]
当 gcd(Y)=1\gcd(Y)=1gcd(Y)=1 时,gcd(x,Y)=1\gcd(x,Y)=1gcd(x,Y)=1 必定满足,所以我们考虑用 gcd(x,Y)=1\gcd(x,Y)=1gcd(x,Y)=1 的所有情况减去 gcd(Y)=1\gcd(Y)=1gcd(Y)=1 的所有情况。
即上式变为 :
∑Y⊆U∑x⊈Y[gcd(x,Y)=1]−∑Y⊆U∑x⊈Y[gcd(Y)=1]\sum_{Y \subseteq U}\sum_{x\not\subseteq Y}[\gcd(x,Y)=1]-\sum_{Y \subseteq U}\sum_{x\not\subseteq Y}[\gcd(Y)=1]Y⊆U∑x⊆Y∑[gcd(x,Y)=1]−Y⊆U∑x⊆Y∑[gcd(Y)=1]
先考虑前半部分,当我们从一个 gcd=1\gcd=1gcd=1 的集合 Z(∣Z∣>1)Z(|Z|>1)Z(∣Z∣>1) 取出任何一个元素 xxx,都满足 gcd(x,Z−x)=1\gcd(x,Z-x)=1gcd(x,Z−x)=1,也就是说 Z−xZ-xZ−x 就是我们要找的集合 YYY。
所以所有的集合产生的贡献就可以表述成:
∑Y⊆U∑x⊈Y[gcd(x,Y)=1]=∑∣Z∣>1,Z⊆U[gcd(Z)=1]∣Z∣\sum_{Y \subseteq U}\sum_{x\not\subseteq Y}[\gcd(x,Y)=1]=\sum_{|Z|>1,Z\subseteq U}[\gcd(Z)=1]|Z|Y⊆U∑x⊆Y∑[gcd(x,Y)=1]=∣Z∣>1,Z⊆U∑[gcd(Z)=1]∣Z∣
对于一个集合 ZZZ,贡献总共 ∣Z∣|Z|∣Z∣ 种构造 xxx 和 YYY 的方法。
我们令 mmm 为 aia_iai 中为 ddd 的倍数的数的数量。
∑∣Z∣>1,Z⊆U[gcd(Z)=1]∣Z∣=∑∣Z∣>1,Z⊆U∑d∣gcd(Z)μ(d)∣Z∣=∑d=1max(ai)μ(d)∑∣Z∣=2mCm∣Z∣∗∣Z∣\begin{aligned}\sum_{|Z|>1,Z\subseteq U}[\gcd(Z)=1]|Z|&=\sum_{|Z|>1,Z\subseteq U}\sum_{d|\gcd(Z)}\mu(d)|Z| \\&=\sum_{d=1}^{max(a_i)}\mu(d)\sum_{|Z|=2}^{m}C_{m}^{|Z|}*|Z| \end{aligned}∣Z∣>1,Z⊆U∑[gcd(Z)=1]∣Z∣=∣Z∣>1,Z⊆U∑d∣gcd(Z)∑μ(d)∣Z∣=d=1∑max(ai)μ(d)∣Z∣=2∑mCm∣Z∣∗∣Z∣
∑∣Z∣=2mCm∣Z∣∗∣Z∣\sum_{|Z|=2}^{m}C_{m}^{|Z|}*|Z|∑∣Z∣=2mCm∣Z∣∗∣Z∣ 可以表述成 ∑i=2nCni∗i\sum_{i=2}^{n}C_{n}^i*i∑i=2nCni∗i,我们令 f(n)=∑i=1nCni∗if(n)=\sum_{i=1}^{n}C_{n}^i*if(n)=∑i=1nCni∗i。
则上述式子写成:
∑d=1max(ai)μ(d)(f(m)−m)\sum_{d=1}^{max(a_i)}\mu(d)(f(m)-m)d=1∑max(ai)μ(d)(f(m)−m)
根据组合数学公式 (1+x)n=∑i=0nCni∗xi(1+x)^n=\sum_{i=0}^nC_n^i*x^i(1+x)n=∑i=0nCni∗xi,两边求导可得 ∑i=0nCni∗i=n∗2n−1\sum_{i=0}^{n}C_{n}^i*i=n*2^{n-1}∑i=0nCni∗i=n∗2n−1 (而 i=0i=0i=0 时无贡献可省去,即为上式)。
则有 f(x)=2f(x−1)+2x−1f(x)=2f(x-1)+2^{x-1}f(x)=2f(x−1)+2x−1
再考虑后半部分,道理类似,任意集合 YYY,产生的贡献为 n−∣Y∣n-|Y|n−∣Y∣,n−∣Y∣n-|Y|n−∣Y∣ 即为 xxx 的数量:
∑Y⊆U∑x⊈Y[gcd(Y)=1]=∑Y⊆U[gcd(Y)=1](n−∣Y∣)=∑Y⊆U∑d∣gcd(Y)μ(d)∗(n−∣Y∣)=∑d=1max(ai)μ(d)((n∗∑∣Y∣=1mCmi)−(∑∣Y∣=1mCmi∗∣Y∣))=∑d=1max(ai)μ(d)(n∗(2m−1)−f(m))\begin{aligned}\sum_{Y \subseteq U}\sum_{x\not\subseteq Y}[\gcd(Y)=1]&=\sum_{Y \subseteq U}[\gcd(Y)=1](n-|Y|)\\&=\sum_{Y \subseteq U}\sum_{d | \gcd(Y)}\mu(d)*(n-|Y|)\\&=\sum_{d=1}^{max(a_i)}\mu(d)((n*\sum_{|Y|=1}^mC_m^i)-(\sum_{|Y|=1}^mC_m^i*|Y|))\\&=\sum_{d=1}^{max(a_i)}\mu(d)(n*(2^m-1)-f(m))\end{aligned}Y⊆U∑x⊆Y∑[gcd(Y)=1]=Y⊆U∑[gcd(Y)=1](n−∣Y∣)=Y⊆U∑d∣gcd(Y)∑μ(d)∗(n−∣Y∣)=d=1∑max(ai)μ(d)((n∗∣Y∣=1∑mCmi)−(∣Y∣=1∑mCmi∗∣Y∣))=d=1∑max(ai)μ(d)(n∗(2m−1)−f(m))
合并起来可得:
∑d=1max(ai)μ(d)(f(m)−m)−∑d=1max(ai)μ(d)(n∗(2m−1)−f(m))=∑d=1max(ai)μ(d)(2f(m)−m−n∗(2m−1))\sum_{d=1}^{max(a_i)}\mu(d)(f(m)-m)-\sum_{d=1}^{max(a_i)}\mu(d)(n*(2^m-1)-f(m))=\sum_{d=1}^{max(a_i)}\mu(d)(2f(m)-m-n*(2^m-1))d=1∑max(ai)μ(d)(f(m)−m)−d=1∑max(ai)μ(d)(n∗(2m−1)−f(m))=d=1∑max(ai)μ(d)(2f(m)−m−n∗(2m−1))
其中 f(1)=1,f(x)=2f(x−1)+2x−1f(1)=1,f(x)=2f(x-1)+2^{x-1}f(1)=1,f(x)=2f(x−1)+2x−1。
这里有一个小 tricktricktrick,我们不需要枚举 aia_iai 的因子来计算贡献。而是通过统计每个 aia_iai 的贡献,然后利用枚举 ddd 的同时,统计它的倍数的贡献(调和级数复杂度 nlognnlognnlogn,如果枚举 aia_iai 的因子则是 nain\sqrt{a_i}nai 复杂度)
#include<bits/stdc++.h>
using namespace std;
#define IOS ios::sync_with_stdio(false); cin.tie(0); cout.tie(0)
#define endl '\n'
typedef long long LL;
const int maxn = 1e7 + 1;
const LL mod = 1e9 + 7;
int bit[maxn], f[maxn], d[maxn], mu[maxn], prime[maxn], cnt;
bool judge[maxn];
inline void init(){
mu[1] = 1;
for(int i = 2; i < maxn; i++){
if(!judge[i]){
prime[++cnt] = i;
mu[i] = -1;
}
for(int j = 1, v; j <= cnt && (v = i * prime[j]) < maxn; j++){
judge[v] = true;
if(!(i % prime[j])) break;
mu[v] = -mu[i];
}
}
f[1] = 1;
bit[0] = 1;
bit[1] = 2;
for(int i = 2; i < maxn; i++){
f[i] = (2LL * f[i - 1] + bit[i - 1]) % mod;
bit[i] = bit[i - 1] * 2LL % mod;
}
}
int main(){
IOS;
init();
int n, a;
cin >> n;
for(int i = 1; i <= n; i++){
cin >> a;
d[a]++;
}
int ans = mu[1] * (2LL * f[n] % mod - n - n * (bit[n] - 1LL) % mod) % mod;
for(int i = 2; i < maxn; i++){
if(mu[i]){
int tol = 0;
for(int j = i; j < maxn; j += i) tol += d[j];
ans = (ans + mu[i] * (2LL * f[tol] % mod - tol - n * (bit[tol] - 1LL) % mod)) % mod;
}
}
cout << (ans + mod) % mod << endl;
}