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 。

这是我做过的最可（du" role="presentation" style="position: relative;">dudu）做（liu" role="presentation" style="position: relative;">liuliu）的NOIP" role="presentation" style="position: relative;">NOIPNOIP模拟级别的数论题。考场上打算写随机算法，想了想30" role="presentation" style="position: relative;">3030暴力更稳，于是交了暴力，结果10" role="presentation" style="position: relative;">1010分滚粗了。

这题让我们联想到正难则反的思想，题目上要我们求出所有可能解的方案数，那么我们这样想，我们先把总方案数求出来，显然是n∗(2n−1−1)" role="presentation" style="position: relative;">n∗(2n−1−1)n∗(2n−1−1)。

然后我们将不合法的情况去掉，我们设加入的数为number" role="presentation" style="position: relative;">numbernumber，那么我们要去掉的就是在加入number" role="presentation" style="position: relative;">numbernumber之前就已经使得gcd" role="presentation" style="position: relative;">gcdgcd为1" role="presentation" style="position: relative;">11的集合对应的number" role="presentation" style="position: relative;">numbernumber的选择方案数之和，以及在加入number" role="presentation" style="position: relative;">numbernumber之后gcd" role="presentation" style="position: relative;">gcdgcd仍然不为1" role="presentation" style="position: relative;">11的方案数。这东西怎么算啊？变形代数式看看吧。

首先解决第二种情况，即在加入number" role="presentation" style="position: relative;">numbernumber之后gcd" role="presentation" style="position: relative;">gcdgcd仍然不为1" role="presentation" style="position: relative;">11的方案数。这种情况直接统计所有集合中使得gcd" role="presentation" style="position: relative;">gcdgcd不为一的方案数即可，即用总方案数减去所有集合中使得gcd" role="presentation" style="position: relative;">gcdgcd为1" role="presentation" style="position: relative;">11的方案数。

另外，我们还知道，情况二其实是让我们求这个东西：∑|S|,gcd(S)!=1|S|" role="presentation" style="position: relative;">∑|S|,gcd(S)!=1|S|∑|S|,gcd(S)!=1|S|

这样的话，我们似乎将第二种情况转化成了第一种情况的子问题。怎么求出所有集合中使得gcd" role="presentation" style="position: relative;">gcdgcd为1" role="presentation" style="position: relative;">11的方案数呢？

考虑容斥原理，我们先将gcd" role="presentation" style="position: relative;">gcdgcd为1" role="presentation" style="position: relative;">11的倍数的集合数求出，减去gcd" role="presentation" style="position: relative;">gcdgcd为2" role="presentation" style="position: relative;">22的倍数的集合数……一直这样迭代下去，然后有个显然的结论，每个gcd" role="presentation" style="position: relative;">gcdgcd的容斥系数就对应着它们的莫比乌斯函数值。于是只需要提前用线性筛筛出所有数的莫比乌斯函数值，每次将情况二的答案加上mu[i]∗(2sum[i]−1)" role="presentation" style="position: relative;">mu[i]∗(2sum[i]−1)mu[i]∗(2sum[i]−1)就行了，其中sum[i]" role="presentation" style="position: relative;">sum[i]sum[i]表示所有数列中所有是i" role="presentation" style="position: relative;">ii的倍数的数的个数。

情况一？考虑用情况二反推情况一，再次将式子变形：我们要求的是在加入number" role="presentation" style="position: relative;">numbernumber之前就已经使得gcd" role="presentation" style="position: relative;">gcdgcd为1" role="presentation" style="position: relative;">11的集合对应的number" role="presentation" style="position: relative;">numbernumber的选择方案数之和，即∑S,gcd(S)==1(n−|S|)" role="presentation" style="position: relative;">∑S,gcd(S)==1(n−|S|)∑S,gcd(S)==1(n−|S|),显然这个式子可以展开成∑S,(gcd)==11−∑S,gcd(S)==1|S|" role="presentation" style="position: relative;">∑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;
}