3月6日 莫比乌斯反演（poj 3904）

莫比乌斯反演可以加速运算，今天大致了解了应用，但是推导还是印象不深的。莫比乌斯反演两种形式，第一种是约数形式：
$f(n)= \sum \limits_{d|n}g(d) \rightarrow g(n)=\sum\limits_{d|n}\mu(d)f(\frac nd)$
第二种是倍数形式：
$f(n)= \sum \limits_{n|d}g(d) \rightarrow g(n)=\sum\limits_{n|d}\mu(\frac dn)f(d)$
其中$g(d)$是我们的目标函数，严格来说这里d是变量。比如说下面说的这道题：
Sky Code
Time Limit: 1000MS Memory Limit: 65536K
Total Submissions: 3017 Accepted: 1021
Description

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.
Input

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.
Output

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

4
2 3 4 5
4
2 4 6 8
7
2 3 4 5 7 6 8
Sample Output

1
0

这里我们构造g(d)就是公因数为d的四元数组组合个数。则f(n)就是公因数为d,2d,3d,…的组合个数之和，这个是比较容易求得的，可以用以下方式筛：对d=i，遍历i，2i，…（i<=max），如果读取的数中有这个数则count[i]++，表示含有因子d的数字个数，从中取出4个数（计算组合数即可）。
那么用第二种形式就可以解决问题了。注意莫比乌斯筛，是一个线性筛。

#include<cstdio>
#include<algorithm>
#define N 10000+3
#define ll long long
int mu[N], pri[N];
bool nopri[N];
int cnt;
void ini()
{
    mu[1] = 1;
    for (int i = 2; i < N; i++)
    {
        if (!nopri[i])
        {
            mu[i] = -1;
            pri[cnt++] = i;
            //printf("pri[%d]=%d\n", cnt - 1, i);
        }
        for (int j = 0; j < cnt&&pri[j] * i < N; j++)
        {
            nopri[pri[j] * i] = 1;
            //printf("%d\n", pri[j]);
            if (i%pri[j] == 0)
            {
                mu[i*pri[j]] = 0;
                break;
            }
            else
            {
                mu[i*pri[j]] = -mu[i];
            }
        }
    }
}
ll com(int m)
{
    return (ll)m*(m - 1)*(m - 2)*(m - 3) / 24;
}
bool data[N];//1表示读取到了这个数
//data在main中memset
int np;
int M;
int count[N];
ll F()
{
    ll res = 0;
    for (int i = 1; i <= M; i++)
    {
        for(int j=i;j<=M;j+=i)
            if (data[j])
            {
                count[i]++;
            }
        if (count[i] >= 4)res += (ll)mu[i]*com(count[i]);
    }
    return res;
}
int main()
{
    ini();
    while (scanf("%d", &np)!= EOF)
    {
        memset(data, 0, sizeof(data));
        memset(count, 0, sizeof(count));
        M = 0;
        for (int i = 0; i < np; i++)
        {
            int tmp;
            scanf("%d", &tmp);
            data[tmp] = 1;
            if (tmp > M)M = tmp;
        }
        if (np < 4)
        {
            printf("0\n");
            continue;
        }
        printf("%lld\n", F());
    }
    return 0;
}