Gcd
题目背景
SOURCE：NOIP2015-SHY-2
题目描述
给出n个正整数，放入数组 a 里。
问有多少组方案，使得我从 n 个数里取出一个子集，这个子集的 gcd 不为 1 ，然后我再从剩下的数中取出一个数，把他放进刚刚取出的子集里，使得 gcd 为 1 。
输出方案数 mod （10^9 + 7）。
输入格式
第一行一个数 n 。
第二行 n 个数，表示 a 数组。
输出格式
输出一个数表示答案。
样例数据 1
输入
3
2 3 2
输出
5
备注
【样例说明】
set ：2 number：3 count：2
set ：3 number：2 count：2
set ：2 2 number：3 count：1
total = 2 + 2 + 1 = 5
【数据范围】
对 30% 的输入数据 ：1≤n≤10；
对 100% 的输入数据 ：1≤n≤500000；2≤a中元素≤10000000 。
这是我做过的最可（dudu）做（liuliu）的NOIPNOIP模拟级别的数论题。考场上打算写随机算法，想了想3030暴力更稳，于是交了暴力，结果1010分滚粗了。
这题让我们联想到正难则反的思想，题目上要我们求出所有可能解的方案数，那么我们这样想，我们先把总方案数求出来，显然是n∗(2n−1−1)n∗(2n−1−1)。
然后我们将不合法的情况去掉，我们设加入的数为numbernumber，那么我们要去掉的就是在加入numbernumber之前就已经使得gcdgcd为11的集合对应的numbernumber的选择方案数之和，以及在加入numbernumber之后gcdgcd仍然不为11的方案数。这东西怎么算啊？变形代数式看看吧。
首先解决第二种情况，即在加入numbernumber之后gcdgcd仍然不为11的方案数。这种情况直接统计所有集合中使得gcdgcd不为一的方案数即可，即用总方案数减去所有集合中使得gcdgcd为11的方案数。
另外，我们还知道，情况二其实是让我们求这个东西：∑|S|,gcd(S)!=1|S|∑|S|,gcd(S)!=1|S|
这样的话，我们似乎将第二种情况转化成了第一种情况的子问题。怎么求出所有集合中使得gcdgcd为11的方案数呢？
考虑容斥原理，我们先将gcdgcd为11的倍数的集合数求出，减去gcdgcd为22的倍数的集合数……一直这样迭代下去，然后有个显然的结论，每个gcdgcd的容斥系数就对应着它们的莫比乌斯函数值。于是只需要提前用线性筛筛出所有数的莫比乌斯函数值，每次将情况二的答案加上mu[i]∗(2sum[i]−1)mu[i]∗(2sum[i]−1)就行了，其中sum[i]sum[i]表示所有数列中所有是ii的倍数的数的个数。
情况一？考虑用情况二反推情况一，再次将式子变形：我们要求的是在加入numbernumber之前就已经使得gcdgcd为11的集合对应的numbernumber的选择方案数之和，即∑S,gcd(S)==1(n−|S|)∑S,gcd(S)==1(n−|S|),显然这个式子可以展开成∑S,(gcd)==11−∑S,gcd(S)==1|S|∑S,(gcd)==11−∑S,gcd(S)==1|S|,哇，后面一坨不是可以跟情况二合并吗？好像还可以抵掉一坨东西。最后发现只需要容斥计算乱搞一下就行了啊。
代码如下：
#include<bits/stdc++.h>
#define inf 10000005
#define mod 1000000007
#define ll long long
using namespace std;
int n,tot=0,mu[inf],prime[inf],cnt[inf],sum[inf];
ll mul[inf],a=0,b=0;
bool pr[inf];
inline ll read(){
ll ans=0;
char ch=getchar();
while(!isdigit(ch))ch=getchar();
while(isdigit(ch))ans=(ans<<3)+(ans<<1)+ch-'0',ch=getchar();
return ans;
}
int main(){
n=read();
for(int i=1;i<=n;++i){
int x=read();
++cnt[x];
}
mul[0]=1;
for(int i=1;i<=n;++i)mul[i]=mul[i-1]<<1,mul[i]%=mod;
mu[1]=1;
for(int i=2;i<=inf;++i){
if(!pr[i])prime[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot&&i*prime[j]<=inf;++j){
int k=i*prime[j];
pr[k]=1;
if(i%prime[j]==0){
mu[k]=0;
break;
}
mu[k]=-mu[i];
}
}
for(int i=1;i<=inf;++i)
for(int j=i;j<=inf;j+=i)sum[i]+=cnt[j];
for(int i=1;i<=inf;++i)
if(mu[i]&&sum[i]){
a=(a+mu[i]*mul[sum[i]]-mu[i])%mod;
b=(b+mu[i]*mul[sum[i]-1]*sum[i])%mod;
}
a=a*n%mod;
b=(b<<1)%mod;
a=(b-a+mod)%mod;
printf("%lld",a);
return 0;
}