2018.06.29 NOIP模拟 Gcd(容斥原理)

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(2n11)n∗(2n−1−1)

然后我们将不合法的情况去掉,我们设加入的数为numbernumber,那么我们要去掉的就是在加入numbernumber之前就已经使得gcdgcd11的集合对应的numbernumber的选择方案数之和,以及在加入numbernumber之后gcdgcd仍然不为11的方案数。这东西怎么算啊?变形代数式看看吧。

首先解决第二种情况,即在加入numbernumber之后gcdgcd仍然不为11的方案数。这种情况直接统计所有集合中使得gcdgcd不为一的方案数即可,即用总方案数减去所有集合中使得gcdgcd11的方案数。
另外,我们还知道,情况二其实是让我们求这个东西:|S|,gcd(S)!=1|S|∑|S|,gcd(S)!=1|S|

这样的话,我们似乎将第二种情况转化成了第一种情况的子问题。怎么求出所有集合中使得gcdgcd11的方案数呢?

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

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

转载于:https://www.cnblogs.com/ldxcaicai/p/9738485.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值