POJ3094 Sky Code(莫比乌斯反演)

POJ3094 Sky Code(莫比乌斯反演)

Sky Code

题意

给你\(n\le 10^5\)个数，这些数\(\le 10^5\)，问这些这些数组成的互不相同的无序四元组(a,b,c,d)使得gcd(a,b,c,d)=1的四元组有多少？

解法

枚举一个约数\(k\)，看看总共有多少个数\(S_k=\{x\}\)满足\(k|x\)。那么可以保证(a,b,c,d)有的一个共同的因子是k，这样的四元组的个数就是

\[F(k)={|S_k|\choose 4} \]

这样算会算重，比如枚举到\(k=4\)再枚举到\(k=2\)，这两者的方案显然有重复，加入有一个四元组满足有一个共同约数是4，那么他们一定也可以满足有一个共同约数是2。我们记\(f(x)=\)最大公因数是\(x\)的四元组的数量。上面的那个大\(F(x)\)就表示有一个公因数(不是最大公因数)是\(x\)的四元组的数量

我们数学模型化这个算重的关系：

\[F(x)=\Sigma_{x|d}f(d) \]

这不就是莫比乌斯反演可以解决的嘛 piece of cake

\[f(x)=\Sigma_{x|d}\mu(d)F(\frac d x) \]

那么把\(f(1)\)求出来就好了

Q:你这样不是O(n^2)吗，你怎么实现可以在正确的复杂度内得到每个数所有的因数？

A:开个桶表示每个\(|S_k|\)，枚举\(i\in [2,\sqrt x]\)，把\(i\)和\(x/i\)都丢在桶里计数。复杂度\(O(n^{1.5})\)注意当\(i=x/i\)的时候只算一次！

//@winlere
#include<iostream>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<algorithm>
#include<vector>
using namespace std;  typedef long long ll;
template < class ccf > inline ccf qr(ccf ret){      ret=0;
      register char c=getchar();
      while(not isdigit(c)) c=getchar();
      while(isdigit(c)) ret=ret*10+c-48,c=getchar();
      return ret;
}

const int maxn=1e4+1;
ll c[maxn][5];
int n;
ll ans;
int buk[maxn];
int cnt[maxn];
bool data[maxn];
int mu[maxn];
bool usd[maxn];
vector < int > ve;

inline void pr(){
      usd[1]=1;mu[1]=0;
      for(register int t=2;t<maxn;++t){
	    if(not usd[t]) ve.push_back(t),mu[t]=-1;
	    for(register int i=0,edd=ve.size();i<edd;++i){
		  register int k=ve[i];
		  if(1ll*k*t>maxn)break;
		  usd[k*t]=1;
		  if(t%k==0) break;
		  mu[k*t]=-mu[t];
	    }
      }
}


int main(){
      c[0][0]=1;
      pr();
      for(register int t=1;t<maxn;++t){
	    c[t][0]=1;
	    for(register int i=1;i<=4;++i)
		  c[t][i]=c[t-1][i-1]+c[t-1][i];
      }
      while(~scanf("%d",&n)){
	    ans=c[n][4];
	    memset(buk,0,sizeof buk);
	    for(register int t=1,data;t<=n;++t){
		  ++buk[data=qr(1)];
		  for(register int i=2;i*i<=data;++i)
			if(data%i==0)
			      if(++buk[i],data/i!=i) ++buk[data/i];
	    }
	    for(register int t=1;t<maxn;++t)
		  if(buk[t]>=4&&mu[t])
			ans+=mu[t]*c[buk[t]][4];
	    printf("%lld\n",ans);
      }
      return 0;
}