poj3904-Sky Code-一起揭露莫比乌斯反演和容斥的一致性

hans774882968

于 2021-05-29 23:44:19 发布

阅读量737

点赞数 2

分类专栏：数论&组合数学

本文链接：https://blog.csdn.net/hans774882968/article/details/117392353

版权

算法数学

数论&组合数学专栏收录该内容

18 篇文章 0 订阅

订阅专栏

莫比乌斯反演

https://vjudge.net/contest/441421#problem/A

原问题可以表示为

$\sum_{1 <= i1<i2<i3<i4<=n} [gcd(a[i1],a[i2],a[i3],a[i4]) == 1]$

根据莫比乌斯反演公式

$[val == 1] = \sum_{d|val}mu[d]$

我们把d移到式子外层，得到

$ans = \sum_{d=1}^{lim}mu[d]*\sum_{1<=i1<i2<i3<i4<=n} [d|gcd(a[i1],a[i2],a[i3],a[i4])]$

其中

$lim=max_{i=1}^{n}(a[i])$

于是设数组里是d的倍数的数共x[d]（这表示x是关于d的函数）个。

$ans=\sum_{d=1}^{lim} mu[d]*C(x[d],4)$

开桶记录每个数的个数，则求x[d]的复杂度是O(lim/d)，于是整个做法的复杂度是nlogn。

容斥

考虑正难则反，原问题可写为

$C(n,4)-\sum_{d>=2}[gcd(a[i1],a[i2],a[i3],a[i4]) == d]$

这个问题依旧没法处理，不妨考虑把gcd转化为倍数关系。我们设集合S[d]为每个数都是d的倍数的4元组集合，则减法部分可以等价写为

$|\bigcup_{d=2}^{n}S[d]|$

我们有一个事实：S[d1]包含S[d2]，如果d1是d2的因数。这意味着S[d2]可以忽略。

我们把d看成若干素因子的乘积，把以上事实用在素因子上，就是：所有含有素因子平方的S[d]都可以忽略，这样只需要考虑lim范围内的素数和它们的组合。事实上，这也是莫比乌斯函数定义的来历。莫比乌斯函数不是从天而降的，它基于集合的包含关系，把有素数平方的d给忽略掉。

考虑到d=6的集合S[6]就是S[2]和S[3]的交集，于是我们获得了能够套容斥原理公式的那些集合：S[p]表示每个数都是p的倍数的4元组集合。再观察原式，我们配凑出，对于有奇数个不同素因子的d，容斥系数取-1，偶数个则取+1。

事实上，这样得到的公式，和莫比乌斯反演解法得到的公式完全一样。

我们宣称，莫比乌斯反演是容斥原理的特例，同属于容斥原理特例的还有二项式反演等。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <vector>
#include <cmath>
#include <queue>
#include <map>
#include <set>
#include <algorithm>
//#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
#define rep(i,a,b) for(int i = (a);i <= (b);++i)
#define re_(i,a,b) for(int i = (a);i < (b);++i)
#define dwn(i,a,b) for(int i = (a);i >= (b);--i)

const int N = 1e4 + 5;

int n,c[N];
bool vis[N];vector<int> primes;int mu[N];

template<typename Type>inline void read(Type &xx){
    Type f = 1;char ch;xx = 0;
    for(ch = getchar();ch < '0' || ch > '9';ch = getchar()) if(ch == '-') f = -1;
    for(;ch >= '0' && ch <= '9';ch = getchar()) xx = xx * 10 + ch - '0';
    xx *= f;
}

void init_primes(int lim){
    vis[1] = true;mu[1] = 1;
    rep(i,2,lim){
        if(!vis[i]) primes.push_back(i),mu[i] = -1;
        re_(j,0,primes.size()){
            int p = primes[j];
            if(i > lim / p) break;
            vis[i * p] = true;
            if(i % p == 0) break;
            mu[i * p] = -mu[i];
        }
    }
}

LL C4(int x){
    return 1LL*x*(x-1)*(x-2)*(x-3)/24;
}

int main(int argc, char** argv) {
    init_primes(N-5);
    while(~scanf("%d",&n)){
        memset(c,0,sizeof c);
        int lim = 0;
        rep(i,1,n){
            int x;read(x);c[x]++;lim = max(lim,x);
        }
        LL ans = 0;
        rep(d,1,lim){
            int v = 0;
            for(int j = d;j <= lim;j += d) v += c[j];
            ans += mu[d] * C4(v);
        }
		printf("%lld\n",ans);
    }
    return 0;
}