容斥原理 HDU-5072

https://vjudge.net/contest/183008#problem/C

解题思路：“如果三个数a, b, c不符合条件，那么一定有一对是互质的，有一对是不互质的。不妨令a, b互质，b, c不互质。于是我们可以枚举b来统计答案。在除了b自己的所有数中，要么与b互质，要么与b不互质。假设n个数中有x个与b不互质的数，那么b对答案的贡献就是(x - 1) * (n - x)。注意这里的求出答案之后要除以2，这是因为如果a, c互质，那么可以交换b, c的位置；如果a, c不互质，那么可以交换a, b的位置，每个答案都被计算了两遍。”

这道题是要求互不相同的n个数中有多少abc这样的组合，满足abc两两互质或者两两不互质。三元关系很伤脑筋，转化为二元关系也许能简化问题（这里脑洞要慢慢打开了）。

由于数据范围不大10^5以内，总组合数C(n,3) longlong不会爆。abc两两互质和两两不互质，就对应着两个互质另两个不互质，这两个集合构成了全集U。不妨把前者称为集合A,后者称为集合B，那么A并B等于U，且A交B为空。U的大小为C(n,3)。

如果a，b，c不符合条件，必然有一对互质，一对不互质，不妨设a，b互质，b，c不互质，于是我们可以枚举b来统计所有的三元组：如果a，c互质那么这样的三元组中b，c可以互换位置；如果a，c不互质，那么a，b可以互换位置。每个答案被算了两遍。

所以只要枚举每个b，统计出k个和它不互质的，那么剩下n-1-k个就是和它互质的，那么三元组就有k*(n-1-k)/2种。

如果会快速统计1~n个数中有多少和b不互质的数，那么题目就差不多就做完了…

对于b不超过10^5，质因子的个数不超过6个（2*3*5*7*11*13 *17>10^5）。找出每个因子的数目。利用容斥原理统计出与b不互质的数的综述。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=100010;
int a[maxn],num[maxn],f[maxn][30],nn[maxn];
int isprime[maxn],prime[10000],cnt=0;
int vis[maxn];
void inti()//指数打表
{
for(int i=2;i<maxn;i++)
{
if(!isprime[i])
{
prime[cnt++]=i;
for(int j=i*2;j<maxn;j+=i)
isprime[j]=1;
}
}
}
//很关键的地方
void fen()//用这种方法求因子不会超时，一个一个分解因子会超时，时间是N根号（N）
{
for(int i=1;i<maxn;i++)
{
for(int j=i;j<maxn;j+=i)
{
if(vis[j])a[i]++;
}
}
}
void fenyin(int n,int k)//分解质因子
{
int c=0;
int m=n;
for(int i=0;i<cnt&&prime[i]*prime[i]<=n;i++)
{

if(m%prime[i]==0)
{
f[k][c++]=prime[i];
while(m%prime[i]==0)
m=m/prime[i];
}
}
if(m>1)
f[k][c++]=m;
nn[k]=c;
}
ll rongchi(int k)//容斥
{
ll sum=0;
ll cn=nn[k];
for(int i=1;i<(1<<cn);i++)
{
int x=0,l=1;
for(int j=0;j<cn;j++)
{
if(i&(1<<j))
x++,l=l*f[k][j];
}
if(x&1)sum+=(a[l]-1);
else sum-=(a[l]-1);
}
return sum;
}
int main()
{

inti();
int t;
scanf("%d",&t);
while(t--)
{
memset(a,0,sizeof(a));
memset(f,0,sizeof(f));
memset(nn,0,sizeof(nn));
memset(vis,0,sizeof(vis));
int n;
scanf("%d",&n);
for(int i=0;i<n;i++)
{
scanf("%d",&num[i]);
vis[num[i]]=1;
fenyin(num[i],i);
}
fen();
long long ans=0;
for(int i=0;i<n;i++)
{
ll s=rongchi(i);
ans+=(ll)s*(ll)(n-s-1);
//cout<<s<<endl;
}
ans/=2;
long long ans1=1;
for(int i=1;i<=3;i++)
ans1=ans1*(n-i+1)/i;
printf("%I64d\n",ans1-ans);
}
}