四元组统计（MSKYCODE - Sky Code）

四元组统计

题目信息

题目链接

Luogu P2714、Luogu SPOJ4191、SPOJ MSKYCODE

题目描述

有 \(n\) 个正整数 \(a _ i\)，你要统计有多少个四元组满足 \(\gcd(a _ i, a _ j, a _ k, a _ l) = 1\)。

输入格式

输入包含多组数据。

对于每组数据：第一行一个正整数 \(n\)，接下来一行 \(n\) 个正整数 \(a _ i\)。

输出格式

若干行，每行对应一个输入数据，表示满足要求的四元组的个数。

样例 #1

样例输入 #1

4
2 3 4 5
4
2 4 6 8
7
2 3 4 5 7 6 8

样例输出 #1

1
0
34

数据范围及约定

对于 \(30\%\) 的数据，\(4 ≤ n ≤ 10\)，且数据组数不超过 \(10\)；

对于 \(100\%\) 的数据，\(4 ≤ n ≤ 10000\)，\(1 ≤ a _ i≤ 10000\)，且数据组数不超过 \(100\)。

思路分析

正难则反。

计算恰好为 \(x\)，这个非常困难，但是求 \(\gcd(a_i,a_j,a_k,a_l)\) 是 \(x\) 的倍数处理却很简单，直接记录因子，是 \(x\) 的倍数的方案数就是 \(vis_x\choose4\)。

定义：\(S_i\) 表示 \(\gcd(a_i,a_j,a_k,a_l)\) 为 \(P_i\) 的倍数的四元组集合。

结果就是求：\(|\displaystyle\bigcup_{i} S_i|\)。但是在 \(10000\) 以内的质数大约有 \(1200\) 个，此时 \(O(2^{cnt})\)，肯定是过不了的。此时有提供两种方法：

方法1

我们直接从一个质数开始，往后面延伸 \(\le 10000\) 的组合。

方法2

其实可以直接 \(1\rightarrow 10000\) 枚举，用莫比乌斯函数暴力求就行，一个的贡献为 \(\mu(i){vis_i\choose4}\)。
时间复杂度：\(O(\max(a_i))\)。此处使用线性筛。

Code

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN = 12000;
int T;
int mu[MAXN],prime[MAXN];
bitset<MAXN> is_prime;
void Euler_sieve(){
	mu[1] = 1;
	for(int i = 2;i<MAXN;i++){
		if(!is_prime[i]){
			prime[++prime[0]] = i;
			mu[i] = -1;
		}
		for(int j = 1;j<=prime[0]&&i*prime[j]<MAXN;j++){
			is_prime[i*prime[j]] = 1;
			if(i%prime[j]==0){
				mu[i*prime[j]] = 0;
				break;
			}else{
				mu[i*prime[j]] = -mu[i];
			}
		}
	}
}
int vis[MAXN];
int n,a[MAXN];
signed main(){
	Euler_sieve();
	while(scanf("%lld",&n)!=EOF){
		memset(vis,0,sizeof(vis));
		for(int i = 1;i<=n;i++){
			scanf("%lld",a+i);
			int tmp = a[i];
			for(int i = 1;i*i<=tmp;i++){
				if(tmp%i==0){
					vis[i]++;
					if(i*i!=tmp) vis[tmp/i]++;
				}
			}
		}
		int ans = 0;
		for(int i = 1;i<MAXN;i++){
			ans += mu[i]*vis[i]*(vis[i]-1)*(vis[i]-2)*(vis[i]-3)/4/3/2/1;
		}
		printf("%lld\n",ans);
	}
	return 0;
}

tag

SPOJ、洛谷、题解
数学、容斥原理、莫比乌斯函数