[日常训练] 太空飞船

【问题描述】

Were非常喜欢太空旅行，但由于他是一个苦逼的程序员所以他买不起自己的太空飞船，所以他决定从Heavy那里偷一艘。
“看！有学妹！”Heavy听了，满心欢喜的跑了出去，Were准备开始下手了。
那么现在只有一个问题：Heavy为他的太空飞船设置了非常牛逼的密码，密码由4个正整数组成，已知它们互不相同，且最大公约数为1，Were显然不可能尝试所有的四元集。幸运的，Were通过某种途径，将数的选择缩小至了一定范围内——一个含有N个数的集合，N<=10000，不过显然，它的四元子集还是太多了，Heavy外出的时间有限，Were必须在极短的时间内破译密码，成败在此一举！快帮他看看究竟有多少种可能的子集。

【输入格式】

输入数据一共两行，第一行一个数为N，第二行一行N个数，给出可能为密码的数字组成的集合，输入数据中所有数不超过10000，保证40%的数据中，N<=50。

【输出格式】

一行一个整数表示存在多少种可能的密码。

【输入输出样例】

样例输入1：
4
2 3 4 5
样例输出1：
1

样例输入2：
4
2 4 6 8

样例输出2：
0

样例输入3：
7
2 3 4 5 7 6 8

样例输出3：
34

【数据范围与约定】

对于40%的数据：n<=50。
对于100%的数据：n<=10000。

【40分】：$O(n^4)$暴力枚举

【60分】：动态规划

设$f[i][j][k]$表示选取到第$i$个数，已选取的数的最大公约数为$j$，选取的个数为$k(0 \le k \le 4)$的总方案数。转移方程如下：

$$f[i][Gcd(a[i], j)][k] += f[i - 1][j][k - 1]$$
则最终答案为 $f[n][1][4]$，其中第一维 $i$可用滚动数组优化。

【100分】：容斥原理

先附上容斥原理介绍
然后我们考虑对问题进行转化：
互质的四元集方案数 = 总方案数 - 不互质的四元集方案数
记$num[i]$表示集合中为 $i$的倍数的元素个数（也就是存在因子$i$的元素个数），则 $Answer = C^4_n - \sum C^4_{num[i]}$
但这样明显存在问题，例如 $num[4]$其实已经包含在 $num[2]$里了， $num[6]$其实已经包含在 $num[3]$里了……所以这里的 $i$必然是一个素数；同时可能存在一个元素包含互不相同的多个素数，那么它在$num[i]$中会被多次统计，因此也要扣除。于是根据容斥原理，我们最终得到的

$$Answer = C^4_n - C^4_{num[2]} - C^4_{num[3]} - C^4_{num[5]}-……$$
$$+ C^4_{num[2 \times 3]} + C^4_{num[2 \times 5]} + ……$$
$$- C^4_{num[2 \times 3 \times 5]} - ……$$
$$= C^4_n + \sum C^4_{num[i]} \times (-1)^{cnt}$$
其中 $i(1 \le i \le 10000)$为一个或多个互不相同的素数的乘积， $cnt$表示 $i$所包含的素数个数。
最后我们利用上述推论来统计答案，总复杂度为$O(n\sqrt n)$。

【代码】

#include <cstdio>
#include <iostream>

using namespace std;
const int N = 1e4, M = N + 5;

int a[M], num[M], n;
long long Ans, c[M];

inline int get()
{
    char ch; int res = 0; bool f = false;
    while (((ch = getchar()) < '0' || ch > '9') && ch != '-');
    if (ch == '-') f = true;
     else res = ch - '0';
    while ((ch = getchar()) >= '0' && ch <= '9')
     res = (res << 3) + (res << 1) + ch - '0';
    return f ? -res : res;
}

int main()
{
    freopen("spacecraft.in", "r", stdin);
    freopen("spacecraft.out", "w", stdout);
    for (int i = 4; i <= N; ++i)
     c[i] = (long long)i * (i - 1) * (i - 2) * (i - 3) / 24;
    //预处理出组合数C 
    n = get();   
    for (int i = 1; i <= n; ++i) num[a[i] = get()] = 1;
    for (int i = 1; i <= n; ++i)
     for (int j = 2; j * j <= a[i]; ++j)
     if (a[i] % j == 0)
     {
        num[j]++;
        if (j * j != a[i]) num[a[i] / j]++;
     }
    Ans = c[n];
    for (int i = 1; i <= N; ++i)
    {
        bool f = false;
        for (int j = 2; j * j <= i; ++j)
         if (i % j == 0 && i / j % j == 0)
         {
            f = true;
            break;
         }
        if (f) continue;
        int cnt = 0, x = i;
        for (int j = 2; j * j <= x; ++j)
         if (x % j == 0) 
         {
            cnt++;
            while (x % j == 0) x /= j;
         }
        if (x != 1) cnt++;
        if (cnt & 1) Ans -= c[num[i]];
         else Ans += c[num[i]];
    }
    cout << Ans << endl;
    fclose(stdin); fclose(stdout);
    return 0;
}