Coprime HDU - 5072【容斥定理+状态压缩】

题目

给出 \(n\) 个互不相同的数，求出满足要么 \(3\) 个数都两两互质，要么都不两两互质的三元组的个数。(\(3\leq n \leq 10^5,1\leq a_i \leq 10^5\))
题目链接：https://vjudge.net/problem/HDU-5072

分析

如果直接从正面考虑怎么求解的话，会发现难以处理。这时，需要逆向思维考虑。可以发现，不满足要求的三元组的一定存在其中两个数互质而另两个数不互质。假设：\(a,b\) 互质，\(b,c\) 不互质。可以通过枚举 \(b\) 来求解。假设与 \(b\) 不互质的数的个数为 \(x\) ，那么 \(b\) 对答案的贡献为 \(\frac{(x-1)(n-x)}{2}\)。之所以要除 \(2\) ，是因为如果 \(a,c\) 互质，那么可以交换 \(b,c\) 的位置；如果 \(a,c\) 不互质，那么可以交换 \(a,b\) 的位置，每个答案都被计算了两遍。

接下来，问题关键在于如何求出有多少个数与 \(b\) 不互质。可以先利用埃式筛筛出每个数的倍数的个数（在所给数中），然后，对枚举的每个 \(b\) 进行质因子分解，得到其各个质因子。利用容斥定理，就可以求出与 \(b\) 不互质的数的个数。最后，要总的个数：\(C(n,3)\) 减去即可。

代码

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int N=1e5+5;
int a[N],prime[N],cnt,num[N],p[10];
bool vis[N],book[N];
void init()
{
    int maxn=1e5;
    for(int i=2;i<=maxn;i++)
    {
        if(!vis[i]) prime[++cnt]=i;
        for(int j=1;j<=cnt&&prime[j]*i<=maxn;j++)
        {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}
void solve()
{
    int maxn=1e5;
    for(int i=1;i<=maxn;i++)
        num[i]=0;
    for(int i=1;i<=maxn;i++)
    {
        for(int j=i;j<=maxn;j+=i)
            if(book[j]) num[i]++;
    }
}
int devide(int x)
{
    int res=0;
    for(int i=1;prime[i]*prime[i]<=x&&i<=cnt;i++)
    {
        if(x%prime[i]==0)
        {
            p[res++]=prime[i];
            while(x%prime[i]==0)
                x/=prime[i];
        }
    }
    if(x>1)
        p[res++]=x;
    return res;
}
int main()
{
    int t,n;
    init();
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        for(int i=1;i<=1e5;i++) book[i]=0;
        for(int i=1;i<=n;i++)
        {
            scanf("%d",&a[i]);
            book[a[i]]=1;
        }
        solve();
        ll ans=0;
        for(int i=1;i<=n;i++)//枚举b
        {
            int m=devide(a[i]);
            ll res=0;//和b不互质的数的个数
            for(int j=1;j<(1<<m);j++)
            {
                int cot=0;
                ll tmp=1;
                for(int k=0;k<m;k++)
                {
                    if(j&(1<<k))
                    {
                        tmp*=p[k];
                        cot++;
                    }
                }
                if(cot&1) res+=num[tmp];
                else res-=num[tmp];
            }
            if(m==0) continue;
            ans+=(res-1)*(n-res);
        }//cout<<"ans="<<ans<<endl;
        ans=1LL*n*(n-1)*(n-2)/6-ans/2;
        printf("%lld\n",ans);
    }
    return 0;
}