【已整理】容斥原理学习
容斥原理的题，经常碰到两种解法，一种是质因数分解后采用dfs进行容斥，另一种是因数分解后递推容斥。两种方法在很多情况下实际上都是可以采用的，只是如果dfs可解，那么dfs的方法会快一些，因为因数分解通常在递推的时候是平方的复杂度，其中需要判断因数之间是否存在倍数关系，这是一个本质上的缺陷。但是质因数分解后dfs就不存在这个问题。但是质因数分解的方法在很多情况下并不适用。
1.hdu6106 Classes
不用管原始题意，给出一个简化的题意：给你七个数，分别表示包含性质a的人有多少个，包含性质b的人有多少个，包含性质c的人有多少个，包含性质ab的人有多少个，包含性质bc的人有多少个，包含性质ac的人有多少个。让你求出至少包含一个性质的人有多少个。
一种方法是直接用容斥原理，a+b+c-ab-bc-ac+abc。但是题目还有一个条件，就是给出的这七个数可能不合法，所以这样算无法判断。然后需要取判断什么情况不合法，一开始想到的就是交集的元素个数必然会小于等于交之前每个集合的元素个数，诸如此类的，但是这样做可能包含的不够全面，最简单最暴力的方法就是，直接求出只包含性质a的人有多少个，只包含性质b的人有多少个，只包含性质c的人有多少个，只包含性质ab的人有多少个，只包含性质bc的人有多少个，只包含性质ac的人有多少个，只包含性质abc的人有多少个，然后判断一下这七个数是不是都大于等于0，如果是则说明合法，答案就是这七个数的和。
那么如何求只包含性质a的人的个数呢？
还是用容斥原理，为a-ab-ac+abc，其它同理
#include<iostream> #include<algorithm> #include<cstring> #include<queue> #include<cstdio> using namespace std; int main(){ int t,n; int A,B,C,AB,BC,AC,ABC; scanf("%d",&t); for(int cas=1;cas<=t;cas++){ int res=0; scanf("%d",&n); for(int i=0;i<n;i++){ scanf("%d%d%d%d%d%d%d",&A,&B,&C,&AB,&BC,&AC,&ABC); int a=A-AB-AC+ABC; int b=B-AB-BC+ABC; int c=C-BC-AC+ABC; int ab=AB-ABC; int bc=BC-ABC; int ac=AC-ABC; int abc=ABC; if(a>=0&&b>=0&&c>=0&&ab>=0&&bc>=0&&ac>=0) res=max(res,a+b+c+ab+bc+ac+abc); } printf("%d\n",res); } return 0; }
2.hdu4135 Co-prime
让你求一个区间内与n互质的数的个数，首先本题可以转化为1-x内与n互质的数
一开始就想到了欧拉函数，但是欧拉函数求的是与1-n范围内与n互质的数，显然与本题不符合。
求1-x内与n互质的数是可以用容斥原理来求的。我们先对n进行素因数分解，获得所有的素因数。
一开始我觉得这样算不对，但是由于pi都是质数，所以这么算是完全正确的。
所以算法分为以下几个步骤
1.预处理33333范围内的质数
2.读入一个n，对n进行素因数分解
3.根据容斥原理采用dfs编写求解函数
4.输出答案
5.如果还有输入，转2，否则结束
下面稍微讲解以下dfs函数
ll ans; void dfs(int u,ll div,int flag,ll num){ if(u==tot){ ans+=(num/div)*flag; return; } dfs(u+1,div,flag,num); dfs(u+1,div*pp[u],flag*(-1),num); } ll coprime(ll num){ if(num<=0) return 0; ans=0; dfs(0,1,1,num); return ans; }
由于容斥原理是加减交替的，且每多一项，加号就变减号，减号就变加号。我们看dfs中的递归调用，如果乘上了pp[u]，则说明多了一项，那么flag就要变号，乘以-1即可。如果没有乘上pp[u]，则调用时参数仍为flag
容斥原理的复杂度一般都是组合数的和，也就是2的幂次，也就是子集的个数。
本题中，全集是n的全部素因子，一般不会很多。最极端的情况是，所有的素因子都只出现一次，那也不会很多，因为前十个素数的积差不多就到int类型的上界了。1e3的复杂度绰绰有余。
#include<iostream> #include<algorithm> #include<cstring> #include<queue> #include<cstdio> using namespace std; typedef long long int ll; const int MAXN=33333; bool vis[MAXN]; int prime[MAXN]; int cnt; int pp[MAXN]; int tot; void get_prime(){ memset(vis,false,sizeof(vis)); cnt=0; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[cnt++]=i; for(int j=i+i;j<MAXN;j+=i) vis[j]=true; } } } void decap(int num){ tot=0; for(int i=0;i<cnt&&prime[i]*prime[i]<=num;i++){ if(num%prime[i]) continue; pp[tot++]=prime[i]; while(num%prime[i]==0) num/=prime[i]; } if(num>1) pp[tot++]=num; } ll ans; void dfs(int u,ll div,int flag,ll num){ if(u==tot){ ans+=(num/div)*flag; return; } dfs(u+1,div,flag,num); dfs(u+1,div*pp[u],flag*(-1),num); } ll coprime(ll num){ if(num<=0) return 0; ans=0; dfs(0,1,1,num); return ans; } int main(){ int t; ll a,b; int n; get_prime(); fun(); scanf("%d",&t); for(int cas=1;cas<=t;cas++){ scanf("%lld%lld%d",&a,&b,&n); decap(n); printf("Case #%d: %lld\n",cas,coprime(b)-coprime(a-1)); } return 0; }
3.hdu 5514 Frogs
所有的gi的情况是取决于输入的，如果我们直接在gi上进行容斥原理，情况会有复杂，编写难度和时间复杂度都会有点高。
那我们考虑利用m的因数，因为最大公因数一定是m的因数，所以我们先对m进行因数分解（注意不是素因数分解）
然后每输入一个ai，我们求出ai和m的最大公因数gi，然后对m的因数中，是gi的倍数的打上标记。
本题的容斥有两种写法，一种是常见的写法，另一种，需要维护一个贡献度数组，我在网上搜到的题解全是后者，一度让我以为前者是不可行的。
我首先介绍一下前者，先放出代码。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=33333; int fac[MAXN<<1]; int cnt; void decap(int num){ cnt=0; for(int i=1;i*i<=num;i++){ if(num%i) continue; fac[cnt++]=i; if(i*i!=num) fac[cnt++]=num/i; } } int gcd(int a,int b){ if(a<b) swap(a,b); if(b==0) return a; return gcd(b,a%b); } ll f[MAXN<<1]; int main(){ int t,n,m; int p; scanf("%d",&t); for(int cas=1;cas<=t;cas++){ memset(f,0,sizeof(f)); scanf("%d%d",&n,&m); decap(m); sort(fac,fac+cnt); for(int i=0;i<n;i++){ scanf("%d",&p); int g=gcd(p,m); for(int i=0;i<cnt;i++){ if(fac[i]%g==0) f[i]=1; } } ll res=0; for(int i=cnt-1;i>=0;i--){ if(!f[i]) continue; int tt=m/fac[i]; f[i]=1LL*fac[i]*(tt-1)*tt/2; for(int j=i+1;j<cnt;j++){ if(fac[j]%fac[i]==0) f[i]-=f[j]; } res+=f[i]; } printf("Case #%d: %lld\n",cas,res); } return 0; }
首先，我求出每个数和m的最大公因数g，然后在m的因数中，所有是g的倍数的打上标记。一开始我认为只有刚好是g的因子才要打上标记，事实证明，这样是不对的，那么为什么我们要把所有是g的倍数的因子都打上标记呢？因为无论我们按照怎样的顺序进行容斥，所有符合条件的数的倍数都要算，所以虽然本质上我们求的是所有g的倍数，但是在容斥的时候，我们也要求出g的倍数的倍数，所以在一开始我们需要对g的倍数打上标记。
然后还有一点，就是我们是算出贡献之后，再去容斥，而不是在统计个数的数组上进行容斥，这一点非常好理解，在有些情况下这个顺序可能无所谓，但是大部分情况下，算出贡献之后再去容斥才是正确的选择。
再放出后者的代码。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=33333; int prime[MAXN]; bool vis[MAXN]; int tot; void init(){ memset(vis,false,sizeof(vis)); tot=0; for(int i=2;i<MAXN;i++){ if(!vis[i]) prime[tot++]=i; for(int j=0;j<tot;j++){ if(1LL*i*prime[j]>=1LL*MAXN) break; vis[i*prime[j]]=true; if(i%prime[j]==0) break; } } } int fac[MAXN<<1]; int cnt; void decap(int num){ cnt=0; for(int i=1;i*i<=num;i++){ if(num%i) continue; fac[cnt++]=i; if(i*i!=num) fac[cnt++]=num/i; } } int a[MAXN]; int gcd(int a,int b){ if(a<b) swap(a,b); if(b==0) return a; return gcd(b,a%b); } int f[MAXN<<1]; int main(){ int t,n,m; scanf("%d",&t); for(int cas=1;cas<=t;cas++){ memset(f,0,sizeof(f)); scanf("%d%d",&n,&m); decap(m); sort(fac,fac+cnt); for(int i=0;i<n;i++){ scanf("%d",a+i); int g=gcd(a[i],m); for(int i=0;i<cnt;i++){ if(fac[i]%g==0) f[i]=1; } } ll res=0; for(int i=0;i<cnt;i++){ if(!f[i]) continue; int tt=m/fac[i]; res+=1LL*fac[i]*(tt-1)*tt*f[i]/2; for(int j=i+1;j<cnt;j++) if(fac[j]%fac[i]==0) f[j]-=f[i]; } printf("Case #%d: %lld\n",cas,res); } return 0; }
前者在容斥的时候是隐式容斥，几乎看不到容斥的过程，本代码维护的f数组是一个容斥贡献度数组，看代码吧，非常好理解。
3.poj1091 跳蚤（容斥原理）
题意：给你n，m。表示一张卡片上有n个数，其中必有一个数为m，其它的数是小于等于m的自然数，问你能够组成多少张卡片，使得跳蚤能够跳到距离它为1的左边的那个位置。其实就是贝祖定理的推论扩展到n个数的情形。那么问题就转化为了，有多少种组合方式，使得这n个数的最大公因数为1。
我们先对m进行质因数分解，对于单独某一个质因数p，其对于不合法的组合的贡献为（m/p）^N。不合法的意思就是gcd>1的情况。这样对于质因数的两两乘积就会重复考虑，这里往下运用容斥原理即可。由于是枚举质因数的全部组合数，我们还是采用dfs编写
#include <iostream> #include <cstdio> #include <cstring> #include <algorithm> using namespace std; typedef long long int ll; const int MAXN=33333; bool vis[MAXN]; int prime[MAXN]; int cnt; int pp[MAXN]; int tot; void get_prime(){ memset(vis,false,sizeof(vis)); cnt=0; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[cnt++]=i; for(int j=i+i;j<MAXN;j+=i) vis[j]=true; } } } void decap(int num){ tot=0; for(int i=0;i<cnt&&prime[i]*prime[i]<=num;i++){ if(num%prime[i]) continue; pp[tot++]=prime[i]; while(num%prime[i]==0) num/=prime[i]; } if(num>1) pp[tot++]=num; } ll quick(ll n,ll k){ ll res=1; while(k){ if(k&1) res*=n; n*=n; k>>=1; } return res; } ll ans; int n,m; void dfs(int u,int div,int flag){ if(u>=tot){ ans+=quick(m/div,n)*flag; return; } dfs(u+1,div,flag); dfs(u+1,div*pp[u],flag*(-1)); } int main(){ get_prime(); scanf("%d%d",&n,&m); decap(m); ans=0; dfs(0,1,1); printf("%lld\n",ans); return 0; }
本题也可以采用因数分解+隐式容斥进行编写（POJ暂时交不上题，下面的代码还没有经过验证）
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=33333; ll quick(ll n,ll k){ ll res=1; while(k){ if(k&1) res*=n; n*=n; k>>=1; } return res; } int fac[MAXN<<1]; int cnt; void decap(int num){ cnt=0; for(int i=1;i*i<=num;i++){ if(num%i) continue; fac[cnt++]=i; if(i*i!=num) fac[cnt++]=num/i; } sort(fac,fac+cnt); } ll f[MAXN<<1]; int main(){ int n,m; while(~scanf("%d%d",&n,&m)){ decap(m); for(int i=cnt-1;i>=0;i--){ f[i]=quick(m/fac[i],n); for(int j=i+1;j<cnt;j++){ if(fac[j]%fac[i]==0) f[i]-=f[j]; } } printf("%lld\n",f[0]); } return 0; }
4.hdu5212 Code
题意：给出n个数，求gcd(a[i], a[j]) * (gcd(a[i], a[j]) - 1)的和（1 <= i, j <= n）
a数组中的数的范围为10000
以下copy了别人的题解
显然，对于一个数 x ，以它为 gcd 的两个数一定都是 x 的倍数。如果 x 的倍数在数列中有 k 个，那么最多有 k^2 对数的 gcd 是 x 。
同样显然的是，对于两个数，如果他们都是 x 的倍数，那么他们的 gcd 一定也是 x 的倍数。
所以，我们求出 x 的倍数在数列中有 k 个，然后就有 k^2 对数满足两个数都是 x 的倍数，这 k^2 对数的 gcd，要么是 x ，要么是 2x, 3x, 4x...
并且，一个数是 x 的倍数的倍数，它就一定是 x 的倍数。所以以 x 的倍数为 gcd 的数对，一定都包含在这 k^2 对数中。
如果我们从大到小枚举 x ，这样计算 x 的贡献时，x 的多倍数就已经计算完了。我们用 f(x) 表示以 x 为 gcd 的数对个数。
那么 f(x) = k^2 - f(2x) - f(3x) - f(4x) ... f(tx) (tx <= 10000, k = Cnt[x])
这样枚举每个 x ，然后枚举每个 x 的倍数，复杂度用调和级数计算，约为 O(n logn)。
说白了，就是枚举每个数作为最大公因数，然后这种计算方法会多算一些东西，这些东西可以用容斥搞定，如果我们从大到小递推，就可以列出一个好像和容斥没有关系的式子，其实还是有隐含的容斥在里面的。为什么这样一来列出来的式子看上去和容斥没有关系呢？也就是我们为什么可以一股脑减掉所有的f（tx）呢？因为我们从大到小递推，所以f（tx）在我用到的时候已经是一个确切的值了！！！
这种类型的容斥原理最重要的是要搞明白容斥的对象到底是什么，我前面说的其实都不对，具体容斥什么还是具体问题具体分析。
本题中，我们统计出了每个因子i的倍数有多少个，比如是x个，那么组合方式就是x平方中，那么容斥的就是组合方式，再用容斥完的组合方式数去乘以最大公因数的贡献。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=10005; const int mod=10007; int cnt[MAXN]; void decap(int num){ for(int i=1;i*i<=num;i++){ if(num%i) continue; cnt[i]++; if(i*i!=num) cnt[num/i]++; } } ll f[MAXN]; int main(){ int n,p; while(~scanf("%d",&n)){ memset(cnt,0,sizeof(cnt)); for(int i=0;i<n;i++){ scanf("%d",&p); decap(p); } ll res=0; for(int i=MAXN-1;i>=1;i--){ if(!cnt[i]) continue; cnt[i]*=cnt[i]; for(int j=i+i;j<MAXN;j+=i) cnt[i]-=cnt[j]; res+=1LL*cnt[i]%mod*i%mod*(i-1)%mod; res%=mod; } printf("%lld\n",res); } return 0; }
别人的方法统计的其实也是每个因子i的倍数的个数。
5.CodeForces - 839D Winter is here
在有了上题的基础之后，再来做这题就显得特别简单了。
我们还是统计出以每个数作为因数的数有多少个，然后用容斥定理去做就特别简单了。
取模的部分不能懈怠，WA了几发。。。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=1000005; const ll mod=1e9+7; ll cnt[MAXN]; ll quick(ll n,ll k){ ll res=1; while(k){ if(k&1) res=(res*n)%mod; n=(n*n)%mod; k>>=1; } return res; } int main(){ int n,p; scanf("%d",&n); memset(cnt,0,sizeof(cnt)); for(int i=0;i<n;i++){ scanf("%d",&p); cnt[p]++; } for(int i=1;i<MAXN;i++){ for(int j=i+i;j<MAXN;j+=i) cnt[i]+=cnt[j]; } ll res=0; for(int i=MAXN-1;i>=2;i--){ if(!cnt[i]) continue; cnt[i]=cnt[i]*quick(2,cnt[i]-1)%mod; for(int j=i+i;j<MAXN;j+=i) cnt[i]=(cnt[i]-cnt[j]+mod)%mod; res+=cnt[i]*i%mod; res%=mod; } printf("%I64d\n",res); return 0; }
6. HYSBZ - 2818 Gcd
给定整数N，求1<=x,y<=N且Gcd(x,y)为素数的数对(x,y)有多少对。
这题仍用前两题的算法也是说得通的，但是会T，题目给的n的范围是1e7，应该是专门为了卡掉这种算法的。
如果仍用这种算法，我们仍旧算出每个数的dp值，最后只统计素数的和即可，但是这样写T了。
再仔细讲一下这种算法吧，这种算法一开始要进行因数分解，统计的从本质意义上来不说，不是每个因子出现的个数，而是这个因子的倍数出现了多少次。
然后枚举每个因数，求出每个因数的倍数的组合数（根据题目不同，计算方法不同），然后减去这个因数的倍数的dp值，就得到了当前这个因数的dp值。
对于本题，最后求所有素数的dp值的和即可。
由于这样的算法会T，所以我们必须考虑用别的方法解决。
搜了一下网上的题解，大致有两种方法，一种使用欧拉函数，一种使用莫比乌斯函数。
先介绍欧拉函数的做法，其实类似的欧拉函数的题我之前做过一题。我们枚举每个素数，求最大公因数为该素数的数有多少对。
我们求出与k互质的数的个数dk，且k*prime[i]<=n，s[i]=sigma(dk*prime[i]),k从1到n/prime[i]，就是最大公因数为该素数的有序对个数，由于题目求的是无序对，所以需要乘以2，再减去重复算的，重复算的只有一种，即（prime[i],prime[i]）。两个相等的数的最大公因数就是它们自身，所以重复算的只有这一种情况，我们对所有小于等于n的素数统计一下s[i]，s[i]的和就是最终的答案。
同时注意和1互质的数有1，这对于一个素数p来说，就是p和p的最大公因数是p，这个无论是有序对还是无序对都只有一种，所以这里需要特别处理一下，即我们将phi[1]设置为0，这样算前缀和的时候就不会算进去，其他的前缀和直接乘以2就是有序对的个数了，对于1的情况，单独考虑，即每次考虑一个素数，贡献加1。
莫比乌斯反演根据细节的不同可以写出下列三份代码。具体反演方法见第8题题解。
①直接根据反演公式算
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=10000005; int prime[MAXN]; bool vis[MAXN]; int tot; int mu[MAXN]; void init(){ tot=0; memset(vis,false,sizeof(vis)); mu[1]=1; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[tot++]=i; mu[i]=-1; } for(int j=0;j<tot&&i*prime[j]<MAXN;j++){ vis[i*prime[j]]=true; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else{ mu[i*prime[j]]=0; break; } } } } ll solve(int k,int N){ ll res=0; for(int i=k;i<=N;i+=k) res+=1LL*mu[i/k]*(N/i)*(N/i); return res; } int main(){ init(); int n; scanf("%d",&n); ll res=0; for(int i=0;i<tot&&prime[i]<=n;i++) res+=solve(prime[i],n); printf("%lld\n",res); return 0; }
②根据反演公式，转化成最大公因数为1的情形（其实只要将范围缩小相应的倍数即可，这个我觉得并不是所有情况都可以这样，本题分析一下是可以这样的）
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=10000005; int prime[MAXN]; bool vis[MAXN]; int tot; int mu[MAXN]; void init(){ tot=0; memset(vis,false,sizeof(vis)); mu[1]=1; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[tot++]=i; mu[i]=-1; } for(int j=0;j<tot&&i*prime[j]<MAXN;j++){ vis[i*prime[j]]=true; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else{ mu[i*prime[j]]=0; break; } } } } ll solve(int N){ ll res=0; for(int i=1;i<=N;i++) res+=1LL*mu[i]*(N/i)*(N/i); return res; } int main(){ init(); int n; scanf("%d",&n); ll res=0; for(int i=0;i<tot&&prime[i]<=n;i++) res+=solve(n/prime[i]); printf("%lld\n",res); return 0; }
③在②基础上的一个整除优化（从AC时间来说没什么太大的提升）
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=10000005; int prime[MAXN]; bool vis[MAXN]; int tot; int mu[MAXN]; int sum[MAXN]; void init(){ tot=0; memset(vis,false,sizeof(vis)); mu[1]=1; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[tot++]=i; mu[i]=-1; } for(int j=0;j<tot&&i*prime[j]<MAXN;j++){ vis[i*prime[j]]=true; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else{ mu[i*prime[j]]=0; break; } } } sum[0]=0; for(int i=1;i<MAXN;i++) sum[i]=sum[i-1]+mu[i]; } ll solve(int N){ ll res=0; int curend=1; int pre=0; for(int i=1;i<=N;i=curend+1){ int tt=N/i; curend=N/tt; res+=1LL*(sum[curend]-sum[pre])*tt*tt; pre=curend; } return res; } int main(){ init(); int n; scanf("%d",&n); ll res=0; for(int i=0;i<tot&&prime[i]<=n;i++) res+=solve(n/prime[i]); printf("%lld\n",res); return 0; }
7.hdu6053 TrickGCD
ai的范围 1e5
题意非常简单，直接看上面的截图吧，基于最经典的gcd容斥问题修改得来。如果直接套用经典方法，很快可以得到一个n方的算法。
枚举gcd，从大到小递推进行隐式容斥。计算每个GCD i的时候，dp[i]=π(a[j]/i)，根据题目的数据量，这样n方的算法过不了。
然后我就觉得肯定不是这么做的，其实是可以优化的。
除以i，想下取整，也就是说以i为循环节的自然数，想下取整除以i都是得到相同的结果，我们利用这一良好的性质进行优化。先记录一下每个数出现了多少次，然后求这个统计数组的前缀和。然后用n*调和级数的复杂度去计算dp[i]即可。
把统计贡献度的部分和容斥的部分的循环范围给弄混了，看了以前的代码才发现错误。。。
当然是所有的数都要考虑进去咯，所以循环是从0开始的，由于涉及上下界的问题，所以代码里的特判有点多。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const ll mod=1e9+7; const int MAXN=100005; ll cnt[MAXN]; ll f[MAXN]; ll quick(ll n,ll k){ ll res=1; while(k){ if(k&1) res=res*n%mod; n=n*n%mod; k>>=1; } return res; } int main(){ int t; int n; int p; scanf("%d",&t); for(int cas=1;cas<=t;cas++){ memset(cnt,0,sizeof(cnt)); scanf("%d",&n); for(int i=0;i<n;i++){ scanf("%d",&p); cnt[p]++; } for(int i=1;i<MAXN;i++) cnt[i]+=cnt[i-1]; ll res=0; for(int i=MAXN-1;i>=2;i--){ f[i]=-1; for(int j=0;j<MAXN;j+=i){ int v=min(MAXN-1,j+i-1); int u=j-1; int dd; if(u==-1) u=0; dd=cnt[v]-cnt[u]; if(dd==0) continue; if(f[i]==-1) f[i]=1; f[i]*=quick(j/i,dd); f[i]%=mod; } if(f[i]==-1){ f[i]=0; continue; } for(int j=i+i;j<MAXN;j+=i) f[i]=(f[i]+mod-f[j])%mod; res=(res+f[i])%mod; } printf("Case #%d: %lld\n",cas,res); } return 0; }
8.接第6题，学习一下莫比乌斯反演
第一种莫比乌斯反演
原式 : G(n)=sigma(F(d)) (其中d|n)
反演公式: F(n)=sigma(U(n/d)*G(d)) 这里U是一个函数，他是每一项 G(d) 的系数，他的定义见下面
（强烈建议关于U的定义这一段可以先跳过，先认为他是G的系数就行了，可以跳到下面红字位置）
(1).U是一个关于整数的函数
(2).U[x] = 1 当且仅当 x能够分解成偶数个不同质数的乘积 （其中1不能被分解，所以1的分解出的质数个数是0，所以U[1]=1）
(3).U[x] = -1 当且仅当 x能够分解成奇数个不同质数的乘积
(4).U[x] = 0 除开(2),(3)的其他情况
看上面关于U的定义可能有点看晕了，通俗一点的说
对于一个 x , 分解因式过后 有 x=(p1^e1)*(p2^e2)...*(pr^er)
如果 ei中(1<=i<=r)有一个数ei大于1 那么 U[x] = 0;
不然的话 U[x] = (-1)^r
依旧来两个例子（我最喜欢举例子了 = =）
U[1]=1;定义中的说明
U[2]=-1; 分解式 2=2;
U[6]=1; 分解式 6=2*3
U[9]=0; 9=3^2; 出现了e>1
U[12]=0; 12=2^2*3;
第二种莫比乌斯反演，与第一种类似，做题时根据列出的函数，判断属于哪一种
原式 : G(n)=sigma(F(d)) (其中n|d,d<=N)
反演公式: F(n)=sigma(U(d/n)*G(d)) (其中n|d,d<=N) 这里U[x]的计算方式和上面的相同！！
举个具体的问题作为例子
给出a，b 其中 (1<=a,b<=10^6)
求满足条件的 x,y 的对数，使得 1<=x<=a,1<=y<=b,且gcd(x,y) == 1。
其中 (2,3) (3,2) 算两对！（即有序对）
用之前容斥原理的方法可以写出nlogn的算法。
用莫比乌斯反演的时间复杂度取决于莫比乌斯函数打表的时间复杂度，因为最后统计答案的时间复杂度是O（n）的，而比较好写的莫比乌斯函数打表是nlogn的，也有O（n）的打表方法
回到本题
我们容斥原理和莫比乌斯反演都做了一下，代码都贴一下
先是莫比乌斯反演
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=1e5+10; const ll mod=1e9+7; int mu[MAXN]; int getmu(){ for(int i=1;i<MAXN;i++){ int target=(i==1?1:0); int delta=target-mu[i]; mu[i]=delta; for(int j=i+i;j<MAXN;j+=i) mu[j]+=delta; } } int main(){ int a,b; getmu(); while(~scanf("%d%d",&a,&b)){ int bound=min(a,b); ll res=0; for(int i=1;i<=bound;i++) res+=mu[i]*((a/i)*(b/i)); printf("%d\n",res); } return 0; }
然后是容斥原理
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=1e5+10; const ll mod=1e9+7; int mu[MAXN]; int getmu(){ for(int i=1;i<MAXN;i++){ int target=(i==1>?1:0); int delta=target-mu[i]; mu[i]=delta; for(int j=i+i;j<=n;j+=i) mu[j]+=delta; } } int main(){ int a,b; getmu(); while(~scanf("%d%d",&a,&b)){ int bound=min(a,b); ll res=0; for(int i=1;i<=bound;i++){ res+=mu[i]*((a/i)*(b/i)); } printf("%d\n",res); } return 0; }
9.hdu1695 GCD
求1-b范围内的x,1-d范围内的y，有多少对（x,y），使得他们的GCD等于给定的k
很明显的莫比乌斯反演，用容斥也可以做的应该。
基础算法就是前一题的弱化版，只需要求出gcd==k的情况，即F（k），但是本题要求给出的是无序对的个数！
我们先求出N=min（b，d），然后x在1-N范围内，y在1-N范围内的个数求出来记为t1，则前一部分无序对个数为（t1+1）/2，然后x在1-b范围内，y在1-d范围内，求出来的个数记为t2，则后一部分无序对的个数为t2-t1，总的无序对的个数即为(t1+1)/2+t2-t1
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=1e5+10; const ll mod=1e9+7; bool vis[MAXN]; int prime[MAXN],mu[MAXN]; int tot; void init(){ memset(prime,0,sizeof(prime)); memset(mu,0,sizeof(mu)); memset(vis,false,sizeof(vis)); mu[1]=1; tot=0; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[tot++] = i; mu[i] = -1; } for(int j=0;j<tot&&i*prime[j]<MAXN;j++){ vis[i*prime[j]]=1; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else{ mu[i*prime[j]]=0; break; } } } } ll cal(int a,int b,int k){ ll t1=0; int N=min(a,b); for(int d=k;d<=N;d+=k) t1+=(ll)mu[d/k]*(N/d)*(N/d); ll t2=0; for(int d=k;d<=N;d+=k) t2+=(ll)mu[d/k]*(a/d)*(b/d); return (t1+1)/2+t2-t1; } int main(){ int a,b,c,d,k; init(); int t; scanf("%d",&t); for(int cas=1;cas<=t;cas++){ scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); printf("Case %d: %lld\n",cas,k==0?0:cal(b,d,k)); } return 0; }
10.hdu3388 Coprime
给你三个数n，m，k，范围均为1e9，问你第k个与n，m互质的数是什么。这里让你求的是与两个数都互质的数，这和求与一个数互质的数是没有区别的，因为我们可以分别对这两个数进行素因数分解，然后将那些至少出现一次的素因数找出来就行了，所以即使给你很多个数，让你求和这些数都互质的数，也是没有任何区别的。然后就是怎么求第k个的问题了。
这里有个点需要明确，第k个可能会很大，那再大，题目没有特别说明，也不会超过long long int的范围吧。怎么求第k个呢？如果我有一个范围，求这个范围内与给定的数的互质的个数，那我用质因子去进行容斥，很快就能算出来。所以我二分答案，找到最小的那个使得互质的数的个数为k的范围，这个范围就是所求的答案。又没有注意到long long int的问题，一直WA。。。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=31623; const ll mod=1e9+7; bool vis[MAXN]; int prime[MAXN]; int tot; int pp[105]; int cnt; void get_prime(){ tot=0; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[tot++]=i; for(int j=i+i;j<MAXN;j+=i) vis[j]=true; } } } void init(int n,int m){ cnt=0; for(int i=0;i<tot&&prime[i]*prime[i]<=n;i++){ if(n%prime[i]) continue; pp[cnt++]=prime[i]; while(n%prime[i]==0) n/=prime[i]; } if(n>1) pp[cnt++]=n; for(int i=0;i<tot&&prime[i]*prime[i]<=m;i++){ if(m%prime[i]) continue; pp[cnt++]=prime[i]; while(m%prime[i]==0) m/=prime[i]; } if(m>1) pp[cnt++]=m; sort(pp,pp+cnt); cnt=unique(pp,pp+cnt)-pp; } ll res; void dfs(int u,ll div,int flag,ll N){ if(u>=cnt){ res+=(N/div)*flag; return; } dfs(u+1,div,flag,N); dfs(u+1,div*pp[u],flag*(-1),N); } ll solve(ll N){ res=0; dfs(0,1,1,N); return res; } int main(){ int t,m,n,k; get_prime(); scanf("%d",&t); for(int cas=1;cas<=t;cas++){ scanf("%d%d%d",&m,&n,&k); init(n,m); ll x=1,y=(1LL<<62); ll res; while(x<=y){ ll mid=x+y>>1; ll tt=solve(mid); if(tt>k) y=mid-1; else if(tt<k) x=mid+1; else{ res=mid; y=mid-1; } } printf("Case %d: %lld\n",cas,res); } return 0; }
11.hdu5072 Coprime 好题！
题意：n个数，选出三个数，如果三个数两两互质，或者三个数两两不互质，则这个三元组符合条件，问你有多少个符合条件的三元组。n范围1e5，数的范围1e5。
分析：第一眼，想分别算这两种情况有多少个，但是三个数两两互质不太好算；三个数两两不互质无法转化为三个数的最大公因数大于1，例如2*3,3*5,2*5。
然后我就考虑了反面情况，即里面有一对数互质，或者有两对数互质的情况。然后我想了一种容斥的方法，写完后调试发现跪了，会有重复的情况，以后写容斥，想到一种容斥方法，一定要拿不太简单的小样例尝试一下看看对不对。
去搜了题解，都说是什么同色三角形，其实对于本题来说，也就是一个很简单的思想。
以后碰到3！一定要想到枚举其中一个的方法。
本题中，我们枚举每个数，找到与其互质的数有k个，则使另外一个数与其互质，另外一个数与其不互质，很容易计算出共有k*（n-k-1）种组合方式，我们枚举完每个数后求和，则到的就是所有不合法的情况了吗？不是的，需要除以2。
为什么？假设三个数分别为A，B，C。不失一般性，假设我们枚举的是A，且B与其互质，C与其不互质，则BC互不互质我们不知道，假设BC互质，则交换AC；假设BC不互质，则交换AB，这样得到的结果我想称之为自对称，对于自对称，真正的不重复的结果当然要除以2了。
那么怎么求和当前的数互质的数有多少个呢？
首先在读入阶段，我们需要对每个数进行因数分解，统计出每个因数的倍数有多少个。然后枚举每个数时，我们对当前数进行素因数分解，然后在素因数上进行dfs容斥，很容易计算，注意由于当前数也属于遍历到的因数的倍数之一，所以在计算贡献的时候，当前因数的倍数的个数要减1！！！
时隔多日再做这题，思路清晰了很多。
大体上思路和上述一致，首先正面考虑，发现不好算，所以我们从反面考虑，由于是三个数，所以我们考虑枚举一个数A，另外两个数B和C则要求B和A互质，C和A不互质，假设n个数，和A互质的数有x个，则和A不互质的数有（n-1-x）个，那么贡献是x*（n-1-x）。那么这么算会重复，到底会重复多少呢？看下面的两个图马上就明白了。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=100005; bool vis[MAXN]; int prime[MAXN]; int tot; void get_prime(){ tot=0; memset(vis,false,sizeof(vis)); for(int i=2;i<MAXN;i++){ if(!vis[i]) prime[tot++]=i; for(int j=0;j<tot;j++){ if(i*prime[j]>=MAXN) break; vis[i*prime[j]]=true; if(i%prime[j]==0) break; } } } int fac[MAXN]; int top; void decap(int num){ top=0; for(int i=0;i<tot&&prime[i]*prime[i]<=num;i++){ if(num%prime[i]) continue; fac[top++]=prime[i]; while(num%prime[i]==0) num/=prime[i]; } if(num>1) fac[top++]=num; } ll cnt[MAXN]; int a[MAXN]; ll ans; void dfs(int now,int div,int flag){ if(now>=top){ ans+=flag*(cnt[div]-1); return; } dfs(now+1,div,flag); dfs(now+1,div*fac[now],flag*(-1)); } int main(){ get_prime(); int t,n; scanf("%d",&t); for(int cas=1;cas<=t;cas++){ memset(cnt,0,sizeof(cnt)); scanf("%d",&n); for(int i=0;i<n;i++){ scanf("%d",a+i); cnt[a[i]]++; } for(int i=1;i<MAXN;i++){ for(int j=i+i;j<MAXN;j+=i) cnt[i]+=cnt[j]; } ll res=0; for(int i=0;i<n;i++){ decap(a[i]); ans=0; dfs(0,1,1); res+=ans*(n-1-ans); } ll tot=1LL*n*(n-1)*(n-2)/6; printf("%lld\n",tot-res/2); } return 0; }
12.hdu2841 Visible Trees
题意：就不讲原本题意了，意思就是1-m，1-n中有多少对有序对的最大公因数是1，之前有一题是1-n，1-n，所以可以用欧拉函数做，这题用欧拉函数我想不到怎么做，所以还是用莫比乌斯反演吧，求F（1）。以后碰到让你求最大公因数等于多少的数对有多少对，一定要想到莫比乌斯反演
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=100005; const ll mod=1e9+7; int mu[MAXN]; int getMu(){ for(int i=1;i<MAXN;i++){ int target=(i==1?1:0); int delta=target-mu[i]; mu[i]=delta; for(int j=i+i;j<MAXN;j+=i) mu[j]+=delta; } } ll solve(ll m,ll n){ ll res=0; ll tt=min(m,n); for(int i=1;i<=tt;i++) res+=mu[i]*(m/i)*(n/i); return res; } int main(){ int t,m,n; getMu(); scanf("%d",&t); for(int cas=1;cas<=t;cas++){ scanf("%d%d",&m,&n); printf("%lld\n",solve(m,n)); } return 0; }
13.HYSBZ-2301 Problem b 重要
题意：对于给出的 n 个询问，每次求有多少个数对 (x,y) ，满足 a ≤ x ≤ b ， c ≤ y ≤ d ，且 gcd(x,y) = k ， gcd(x,y) 函数为 x 和 y 的最大公约数。
一开始觉得就是直接用莫比乌斯反演算一下就好了，但是当k比较小时，比如1时，整个算法的时间复杂度是50000的平方的，不能过。所以这种题目是肯定要用到n/i的分块优化的，但是我直接算肯定用不出来分块优化。所以我在这里先用一个小容斥
（1-b，1-d）-（1-a，1-d）-（1-b，1-c）+（1-a，1-c） 区间的开闭问题需要注意一下
然后gcd（x，y）=k也可以化成gcd=1的情况，只要数的范围向下取整除以k即可
莫比乌斯反演最后统计时的分块优化我在上面的题目里也用过，只是这题同时有n和m的限制，所以curend赋值为(n/i)和(m/i)的较小者
#pragma comment(linker, "/STACK:102400000,102400000") #include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=50005; const ll mod=1e9+7; bool vis[MAXN]; int prime[MAXN],mu[MAXN]; int sum[MAXN]; int tot; void init(){ memset(prime,0,sizeof(prime)); memset(mu,0,sizeof(mu)); memset(vis,false,sizeof(vis)); mu[1]=1; sum[0]=0; sum[1]=1; tot=0; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[tot++] = i; mu[i] = -1; } for(int j=0;j<tot&&i*prime[j]<MAXN;j++){ vis[i*prime[j]]=1; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else{ mu[i*prime[j]]=0; break; } } sum[i]=sum[i-1]+mu[i]; } } ll solve(ll n,ll m){ ll res=0; ll tt=min(n,m); ll curend; for(ll i=1;i<=tt;i=curend+1){ curend=min(n/(n/i),m/(m/i)); res+=(sum[curend]-sum[i-1])*(n/i)*(m/i); } return res; } int main(){ ll t,a,b,c,d,k; init(); scanf("%lld",&t); for(int cas=1;cas<=t;cas++){ scanf("%lld%lld%lld%lld%lld",&a,&b,&c,&d,&k); printf("%lld\n",solve(b/k,d/k)+solve((a-1)/k,(c-1)/k)-solve((a-1)/k,d/k)-solve((c-1)/k,b/k)); } return 0; }
14.hdu6134 Battlestation Operational
这题想好好写个完整的题解，所以挂个博文链接
http://www.cnblogs.com/milesgo/articles/7385524.html
时隔两个月又做了一遍，最后询问的处理上又卡了。。。没必要全部预处理出来的啊。。。
我一直在想怎么才能够在规定时间内预处理出来。。。其实询问只有4000个，每次根号n的算一下，就能够时间了。
#include <bits/stdc++.h> using namespace std; typedef long long int ll; const int MAXN=1000005; const ll mod=1e9+7; int prime[MAXN]; bool vis[MAXN]; int tot; int mu[MAXN]; int sum[MAXN]; void mu_init(){ tot=0; memset(vis,false,sizeof(vis)); mu[1]=1; for(int i=2;i<MAXN;i++){ if(!vis[i]){ prime[tot++]=i; mu[i]=-1; } for(int j=0;j<tot&&i*prime[j]<MAXN;j++){ vis[i*prime[j]]=true; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else{ mu[i*prime[j]]=0; break; } } } sum[0]=0; for(int i=1;i<MAXN;i++) sum[i]=(sum[i-1]+mu[i]+mod)%mod; } ll cnt[MAXN]; ll f[MAXN]; ll g[MAXN]; void init(){ mu_init(); //预处理每个数的因子个数 memset(cnt,0,sizeof(cnt)); for(int i=1;i<MAXN;i++){ cnt[i]++; cnt[i]%=mod; for(int j=i+i;j<MAXN;j+=i) cnt[j]++; } f[0]=0; for(int i=1;i<MAXN;i++) f[i]=(f[i-1]+cnt[i-1]+1)%mod; for(int i=1;i<MAXN;i++) f[i]=(f[i]+f[i-1])%mod; } int main(){ int n; init(); while(~scanf("%d",&n)){ ll res=0; int curend; for(int i=1;i<=n;i=curend+1){ int tt=n/i; curend=n/tt; res+=(sum[curend]-sum[i-1]+mod)*f[tt]; res%=mod; } printf("%lld\n",res); } return 0; }
浙公网安备 33010602011771号