POJ 3904 Sky Code 莫比乌斯反演容斥原理

最新推荐文章于 2019-07-27 19:48:21 发布

原创最新推荐文章于 2019-07-27 19:48:21 发布 · 2.3k 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#poj #数学 #acm #算法 #莫比乌斯反演

ACM-解题报告同时被 2 个专栏收录

38 篇文章

订阅专栏

|--数学

22 篇文章

订阅专栏

本文深入探讨了一个数论问题，即在一组给定的整数中，寻找四数的组合，使得这四个数的最大公约数为1。通过容斥原理和莫比乌斯反演的方法，作者提供了求解此类问题的思路，并给出了相应的代码实现。此方法不仅适用于理解数论的基本概念，也展示了算法解决复杂数学问题的灵活性。

原题见POJ 3904

给n个数，求其中四个数的gcd是1的情况有多少种。

从反面考虑，算出gcd不是1的情况，总数取反即是结果。这是容斥原理的思想。当时在做POJ 1091的时候即是这样的想法。画一个vene图，每个集合表示最大公约数为k的倍数情况数。
vene
当k含有素数因子的平方项，如4，12，其实已经被2的情况数覆盖，不必再进行任何处理。只需考虑k是素数的一次方的乘积的情况。当素数个数为奇数个，如2，3，5，30，情况数是向上累加的；而当素数个数是6，10，15时，情况数应该被减去，因为被重复算进过。
其实和莫比乌斯反演里的系数的产生是一样的原理。

有了莫比乌斯反演的知识背景以后，可以这样想：

设 $f(x)$ =四个数的gcd是x的情况数
$F(x)=\sum_{x\mid d}f(d)$ ，即四个数的gcd是x的倍数的情况数
利用上一篇博客提到的莫比乌斯反演的第二种形式

现在要求 $f(1)=\sum_{x\mid d}\mu({d\over x})F( d)|_{x=1}=\sum_d\mu(d)F(d)$ 其中d是1~10000的所有数。

附code

/*--------------------------------------------
 * File Name: POJ 3904
 * Author: Danliwoo
 * Mail: Danliwoo@outlook.com
 * Created Time: 2015-10-06 12:36:28
--------------------------------------------*/

#include <cstdio>
#include <iostream>
#include <cstring>
#include <queue>
#include <algorithm>
#include <cmath>
using namespace std;
#define N 10001
#define LL long long
int mobi[N]={0,1}, num[N], prim[N], cnt = 0, a[N];
void set()
{
    for(int i = 0;i < N;i++) num[i] = 1;
    for(int i = 2;i < N;i++)
    {
        if(num[i])
        {
            prim[cnt++] = i;
            mobi[i] = -1;
        }
        for(int j = 0;j < cnt && prim[j]*i < N;j++)
        {
            num[prim[j]*i] = 0;
            if(i % prim[j] == 0)
            {
                mobi[prim[j]*i] = 0;
                break;
            }
            mobi[prim[j]*i] = -mobi[i];
        }
    }
}
void sep(int n)
{
    int s = sqrt(n);
    for(int i = 1;i <= s;i++) if(n % i == 0)
    {
        a[i]++;
        a[n/i]++;
    }
    if(s*s == n)
        a[s]--;
}
LL get(LL i)
{
    return i*(i-1)*(i-2)*(i-3)/24;
}
int main()
{
    set();
    int n, tp;
    while(~scanf("%d", &n))
    {
        memset(a, 0, sizeof(a));
        for(int i = 0;i < n;i++)
        {
            scanf("%d", &tp);
            sep(tp);
        }
        LL ans = 0;
        for(LL i = 1;i < N;i++)
        {
            if(mobi[i] == 0 || a[i] < 4) continue;
            ans += get(a[i])*mobi[i];
        }
        printf("%lld\n", ans);
    }
    return 0;
}