[题解] SP4191 Sky Code

[题解] SP4191 Sky Code 莫比乌斯反演

题目链接
这是我做的最水的一道莫比乌斯板题。。。不过这题是我第一次真正用到莫比乌斯反演，之前都是用的莫比乌斯函数的性质。。。
我们这里直接定义$f(n)$为$a_1$到$a_n$中，$\gcd (a,b,c,d)=n$的无序四元组个数，$g(n)$为$a_1$到$a_n$中，$\gcd (a,b,c,d)=n的倍数$的无序四元组个数。那么显然有
$$g(n)=\sum _{n|d}f(d)$$
套式子可得
$$f(n)=\sum _{n|d}\mu (d)g\left(\frac{d}{n}\right)$$
我们要求$$f(1)=\sum _{d=1}^{\max a_i}\mu (d)g(d)$$
这个$g$很好求，因为$\gcd (a,b,c,d)=n的倍数$等价于
$$n|a,n|b,n|c,n|d$$我们设$k(i)$为数组中$i$的倍数的个数。
那么容易求得
$$g(n)=\left(\genfrac{}{}{0ex}{}{k(n)}{4}\right)$$
故
$$f(1)=\sum _{d=1}^{\max a_i}\mu (d)\left(\genfrac{}{}{0ex}{}{k(d)}{4}\right)$$
对于$k(d)$，我们用试除法统计每个$a_i$的因子。

#include<bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define pii pair<int,int>
using namespace std;

const double eps = 1e-10;
const double pi = acos(-1.0);
const int maxn = 1e5 + 10;

ll mu[maxn],p[maxn];
ll k[maxn];
ll C4[maxn];
int cnt;
bool isp[maxn];

void getp(){
	isp[1] = mu[1] = 1;
	for(int i = 2; i < maxn; i++){
		if(!isp[i]) p[++cnt] = i,mu[i] = -1;
		for(int j = 1; j <= cnt && i * p[j] < maxn; j++){
			isp[i*p[j]] = 1;
			if(i % p[j]) mu[i*p[j]] = -mu[i];
			else break;
		}
	}
}

void getc(){
	for(ll i = 4; i < maxn; i++) C4[i] = i*(i-1)*(i-2)*(i-3)/24;
}

void cal(int x){
	for(int i = 1; i * i <= x; i++){
		if(x % i == 0){
			k[i]++;
			k[x/i]++;
			if(i * i == x) k[i]--;
		}
	}
}

int n,a[maxn];

void solve(){
	getp();
	getc();
	while(scanf("%d",&n) != EOF){
		memset(a,0,sizeof a);
		memset(k,0,sizeof k);
		int up = 0;
		for(int i = 1; i <= n; i++){
			scanf("%d",&a[i]);
			cal(a[i]);
			up = max(up, a[i]);
		}
		ll ans = 0;
		for(int i = 1; i <= up; i++) ans += mu[i]*C4[k[i]];
		printf("%lld\n",ans);
	}
}

int main()
{
	solve();
	return 0;
}