POJ 3904 Sky Code (容斥＋莫比烏斯反演)

Stancu likes space travels but he is a poor software developer and will never be able to buy his own spacecraft. That is why he is preparing to steal the spacecraft of Petru. There is only one problem – Petru has locked the spacecraft with a sophisticated cryptosystem based on the ID numbers of the stars from the Milky Way Galaxy. For breaking the system Stancu has to check each subset of four stars such that the only common divisor of their numbers is 1. Nasty, isn’t it? Fortunately, Stancu has succeeded to limit the number of the interesting stars to N but, any way, the possible subsets of four stars can be too many. Help him to find their number and to decide if there is a chance to break the system.

In the input file several test cases are given. For each test case on the first line the number N of interesting stars is given (1 ≤ N ≤ 10000). The second line of the test case contains the list of ID numbers of the interesting stars, separated by spaces. Each ID is a positive integer which is no greater than 10000. The input data terminate with the end of file.

For each test case the program should print one line with the number of subsets with the asked property.

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

1 
0 
34

Southeastern European Regional Programming Contest 2008

題目連結：http://poj.org/problem?id=3904

題目大意：給n個不相同的數，從中任意選出4個，使得它們的最大公約數為1，問有多少種選法

題目分析：首先n小於4肯定是0了，C(n,4) = n * (n - 1) * (n - 2) * (n - 3) / 24，本題資料範圍不大，long long即可，差不多是裸的莫比烏斯反演題了，和NOJ 2079幾乎一樣，ac程式碼在POJ上rank 7

#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
int const MAX = 1e4 + 5;
int mob[MAX], p[MAX], cnt[MAX], num[MAX];
bool prime[MAX];
int n, ma;

void Mobius()
{
    memset(prime, true, sizeof(prime));
    int pnum = 0;
    mob[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        if(prime[i])
        {
            p[pnum ++] = i;
            mob[i] = -1;
        }
        for(int j = 0; j < pnum && i * p[j] < MAX; j++)
        {
            prime[i * p[j]] = false;
            if(i % p[j] == 0)
            {
                mob[i * p[j]] = 0;
                break;
            }
            mob[i * p[j]] = -mob[i];
        }
    }
}

ll cal()
{
    for(int i = 1; i <= ma; i++)
        for(int j = i; j <= ma; j += i)
            num[i] += cnt[j];
    ll ans = 0;
    for(int i = 1; i <= ma; i++)
    {
        int x = num[i];
        if(x >= 4)
            ans += (ll) mob[i] * x * (x - 1) * (x - 2) * (x - 3) / 24;
    }
    return ans;
}

int main()
{
    Mobius();
    while(scanf("%d", &n) != EOF)
    {
        ma = 0;
        memset(cnt, 0, sizeof(cnt));
        memset(num, 0, sizeof(num));
        for(int i = 0; i < n; i++)
        {
            int tmp;
            scanf("%d", &tmp);
            cnt[tmp] ++;
            ma = max(ma, tmp);
        }
        if(n < 4)
        {
            printf("0\n");
            continue;
        }
        printf("%lld\n", cal());
    }
}