题目来源：POJ 3904 Sky Code

题意：选出最大公约数为1的四元组的方案

思路：容斥原理 总的方案C(n,4)减去t(1)+t(2)-t(3)+...+(-)^kt(k)

t(i)表示四元组质因子的个数为i的方案数

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxn = 10010;
typedef long long LL;
int a[maxn];
int b[maxn];
int cnt[maxn];
int n, m;
//返回a^p mod n 快速幂
bool vis[maxn];
int prime[maxn]; 

LL C(LL n, LL m)
{
	if(n < m)
		return 0;
	LL ans = 1;
	for(int i = 1; i <= m; i++)
	{
		ans *= n--;
		ans /= i;
	}
	return ans;
}
int main()
{
	while(scanf("%d", &n) != EOF)
	{
		memset(a, 0, sizeof(a));
		for(int i = 0; i < n; i++)
		{
			scanf("%d", &b[i]);
			int x = b[i];
			int sum = 0;
			for(int j = 2; j*j <= x; j++)
			{	
				if(x%j == 0)
				{
					prime[sum++] = j;
					x /= j;
					while(x%j == 0)
						x /= j;
				}
			}
			if(x > 1)
				prime[sum++] = x;
			for(int j = 1; j < (1<<sum); j++)
			{
				int temp = 1;
				int num = 0;
				for(int k = 0; k < sum; k++)
				{
					if(j&(1<<k))
					{
						temp *= prime[k];
						num++;
					}
				}
				a[temp]++;
				cnt[temp] = num;
			}
		}
		LL ans = 0;
		for(int i = 0; i <= 10000; i++)
		{
			if(!a[i])
				continue;
			if(cnt[i]&1)
				ans += C(a[i], 4);
			else
				ans -= C(a[i], 4);
		}
		printf("%I64d\n", C(n, 4)-ans);
	}
	return 0;
}