POJ3904 & SPOJ & 洛谷SP4191 Sky Code 莫比乌斯反演

题目链接（略多）：
SPOJ
POJ3904
洛谷SP4191
vjudge
cn.vjudge.net

题目大意：

给出$n$个数$a_1,a_2,...,a_n$，要求从$n$个数中选出$4$个数$a,b,c,d$，要求$gcd(a,b,c,d)=1.$
求共有多少种取法。
数据范围：$n<=10^4,1<=a_i<=10^4$，每个测试点包含多组数据。

暴力：枚举四个选到的数，然后大力求$gcd...$
时间复杂度$O(n^{4}logn)...$在$n=10^4$的情况下肯定要凉qwq
于是我们珂以考虑上网搜题解……
莫比乌斯反演？
于是我学了一下……（这道题）好像也不是很难？qwq
莫比乌斯反演的主要内容就一条：
若$F(n)=\Sigma {}_{d|n}f(d)$，则有$f(n)=\Sigma {}_{d|n}\mu (d)F(\frac{n}{d})$ ——（1）
但有时候用的是另外一种形式：
若$F(n)=\Sigma {}_{n|d}f(d)$，则有$f(n)=\Sigma {}_{n|d}\mu (\frac{d}{n})F(d)$ ——（2）
如果这里给出详细证明的话会变得很长……所以这里就不证明了……可以看我的博客qwq
但是这两个玄学式子对这道题有什么用呢？
不妨让$f(n)$表示在所有数中选出$4$个，且$4$个数的$gcd$为$n$的选法。
让$F(n)$表示在所有数中选出$4$个，且$4$个数的$gcd$为$n$的倍数的选法。
现在需要求的是$f(1)$。
珂以发现：
$F(n)=\Sigma {}_{n|d}f(d)$
要注意的是，这里$d\le max(a_1,a_2,...,a_n)$，否则$f(d)=0$（最大公因数不可能大于所有数的最大值）
然后把这个式子往（2）一套，得到
$f(n)=\Sigma {}_{n|d}\mu (\frac{d}{n})F(d)$
所以
令$maxa=max(a_1,a_2,...,a_n)$
$f(1)=\Sigma {}_{1|d}\mu (\frac{d}{1})F(d)=\Sigma {}_{d=1}^{maxa}\mu (d)F(d)$
所以如果知道了所有的$\mu (i)和F(i)$，就珂以在$O(n)$的时间复杂度内求解qwq。
然后$\mu (i)$是珂以线性筛的……方法和筛质数差不多……具体见代码qwq
$F(i)$的求法也比较简单，把这些数分解质因数，令$num[i]$表示有多少个数能整除$i$，在分解质因数时求出所有的$num.$
显然，$F(i)=C_{num[i]}^4$
所以$F(i)$能在$O(n\sqrt{n})$复杂度内筛出来。
然后算出答案就珂以了……qwq

代码

#include<stdio.h>
#include<cstring>
#include<algorithm>
#include<math.h>
#include<vector>
#define re register int
#define rl register ll
#define lowbit(x) x&(-x)
using namespace std;
typedef long long ll;
int read() {
	re x=0,f=1;
	char ch=getchar();
	while(ch<'0' || ch>'9') {
		if(ch=='-')	f=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9') {
		x=10*x+ch-'0';
		ch=getchar();
	}
	return x*f;
}
inline void write(int x) {
	if(x>9)	write(x/10);
	putchar(x%10+'0');
}
namespace I_Love {

const int Size=10005;
int n,tot,a[Size],p[Size],prime[Size];
bool vis[Size];
void getp(int maxn) {
	//线性筛莫比乌斯函数
	p[1]=1;
	for(re i=2; i<=maxn; i++) {
		if(!vis[i]) {
			prime[++tot]=i;
			p[i]=-1;
		}
		for(re j=1; j<=tot; j++) {
			int now=i*prime[j];
			if(now>maxn)	break;
			vis[now]=true;
			//如果i能整除prime[j]，就break 
			//可以证明这样i*prime[j+1]也一定会被筛到 
			if(!(i%prime[j]))	break;
			p[now]=-p[i];
		}
	}
}
ll num[Size],F[Size];
void Divide(int x) {
	int sx=sqrt(x);
	for(re i=1; i<=sx; i++) {
		if(!(x%i)) {
			num[i]++;
			num[x/i]++;
		}
	}
	if(sx*sx==x) {
		num[sx]--;
	}
}
void Fujibayashi_Ryou() {
	getp(10000);
	while(scanf("%d",&n)==1) {
		int maxn=0;
		for(re i=1; i<=n; i++) {
			a[i]=read();
			if(a[i]>maxn)	maxn=a[i];
		}
		memset(num,0,sizeof(num));
		for(re i=1; i<=n; i++) {
			//分解质因数，同时求出所有的num
			Divide(a[i]);
		}
		ll ans=0;
		for(re i=1; i<=maxn; i++) {
			F[i]=num[i]*(num[i]-1)*(num[i]-2)*(num[i]-3)/24;
			ans+=F[i]*p[i];
		}
		printf("%lld\n",ans);
	}
}

}
int main() {
	I_Love::Fujibayashi_Ryou();
	return 0;
}