蓝书（算法竞赛进阶指南）刷题记录——POJ3904 Sky Code（容斥原理）

题目：POJ3904.
题目大意：给定$n$个数$a_i$，求有多少个四元组$(w,x,y,z)$满足$gcd(a_w,a_x,a_y,a_z)=1$.
$1\le n,a_i\le 10^4$.

我容斥果然还是太菜了，这种题都要看题解…

首先$gcd(a,b,c,d)=1$不代表两两互质.

考虑把问题转化为用所有四元组的数量$C_n^4$减掉$gcd$不为$1$的四元组数量.

考虑如何求出后者，我们统计$cnt[i]$表示有多少个数是$i$的倍数，那么有$C_{cnt[i]}^4$个四元组的$gcd$是$i$的倍数.

想到这一步就开始考虑容斥，很容易发现$i$的容斥系数即为Mobius函数$\mu (i)$，线性筛预处理即可，时间复杂度$O(n\log n)$.

代码如下：

#include<iostream>
#include<cstdio>
#include<algorithm>
  using namespace std;
  
#define Abigail inline void
typedef long long LL;

const int N=10000;

LL C(LL n){return n<4?0:n*(n-1)*(n-2)*(n-3)/24;}

int pr[N+9],tp,b[N+9];
LL mu[N+9];

void sieve(int n){
  for (int i=2;i<=n;++i) b[i]=1;
  mu[1]=1;
  for (int i=2;i<=n;++i){
  	if (b[i]) pr[++tp]=i,mu[i]=-1;
  	for (int j=1;j<=tp&&i*pr[j]<=n;++j){      //把n写成tp调了好久... 
	  b[i*pr[j]]=0;
	  if (i%pr[j]==0){mu[i*pr[j]]=0;break;}
	  mu[i*pr[j]]=-mu[i];
	}
  }
}

int n,mx;
LL c[N+9],cnt[N+9],ans;

Abigail start(){
  sieve(N);
}

Abigail into(){
  int a;
  for (int i=1;i<=n;++i){
    scanf("%d",&a);
    mx=max(mx,a);
    ++cnt[a];
  }
}

Abigail work(){
  for (int i=1;i<=mx;++i)
    for (int j=i;j<=mx;j+=i)
      c[i]+=cnt[j];
  for (int i=1;i<=mx;++i)
    ans+=mu[i]*C(c[i]);
}

Abigail outo(){
  printf("%lld\n",ans);
  for (int i=1;i<=mx;++i)
    cnt[i]=c[i]=0;
  mx=ans=0;
}

int main(){
  start();
  while (~scanf("%d",&n)){
    into();
    work();
    outo();
  }
  return 0;
}