容斥原理与Mobius函数

·# 容斥原理
百度百科：容斥原理

设$S_1,S_2,\cdots ,S_n$为有限集合，$|S|$表示集合$S$的大小，则：
$$\left|\bigcup _{i=1}^n{S}_i\right|=\sum _{i=1}^n\left|S_i\right|-\sum _{1\le i<j\le n}\left|S_i\cap S_j\right|+\sum _{i\le i<j<k\le k}\left|S_i\cap S_j\cap S_k\right|+\cdots +(-1{)}^{n+1}\left|S_1\cap \cdots \cap S_n\right|$$

多重集

多重集是指包含重复元素的广义集合。设 S = { n 1 ∗ a 1 , n 2 ∗ a 2 , ⋯   , n k ∗ a k } S = \left \{ n_1 * a_1 , n_2 * a_2, \cdots , n_k * a_k \right \} S={n1​∗a1​,n2​∗a2​,⋯,nk​∗ak​}，是由$n_1$个$a_1$，$n_2$个$a_2$$\cdots$$n_k$个$a_k$组成的多重集。则$S$ 的全排列个数为：
$$\frac{n!}{n_1!n_2!\cdots n_k!}$$
设$r\le n_i(\forall \in [1,k]$。从$S$中取出$r$个元素组成一个多重集（不考虑顺序），产生的不同多重集数量为：
$$C_{k+r-1}^{k-1}$$
设$n={\sum }_{i=1}^k{n}_i$，对于任意整数：$r\le n$，从$S$中取出$r$个元素组成的一个多重集（不考虑顺序），产生的不同多重集的数量为：
$$C_{k+r-1}^{k-1}-\sum _{i=1}^k{C}_{k+r-n_i-2}^{k-1}+\sum _{1\le i<j\le k}{C}_{k+r-n_i-n_j-3}^{k-1}-\cdots +(-1{)}^k{C}_{k+r-{\sum }_{i=1}^k{n}_i-(k+1)}$$

例题1：

题目链接

题意：

Devu 有 $N$ 个盒子，第 $i$ 个盒子中有 $A_i$ 枝花。
同一个盒子内的花颜色相同，不同盒子内的花颜色不同。
Devu 要从这些盒子中选出 $M$ 枝花组成一束，求共有多少种方案。
若两束花每种颜色的花的数量都相同，则认为这两束花是相同的方案。
结果需对 $10^9+7$ 取模之后方可输出。

Sol：

根据多重集组合数结论，本体答案为上式：$$C_{k+r-1}^{k-1}-\sum _{i=1}^k{C}_{k+r-n_i-2}^{k-1}+\sum _{1\le i<j\le k}{C}_{k+r-n_i-n_j-3}^{k-1}-\cdots +(-1{)}^k{C}_{k+r-{\sum }_{i=1}^k{n}_i-(k+1)}$$

• tips: 二进制表示集合，利用$lucas$ 求组合数

Code：

#include <bits/stdc++.h>
#define int long long
using namespace std;

const int N = 22, mod = 1e9 + 7;
int a[22], inv[22];
int n ,m;

int qpow(int a, int b)
{
    int res = 1;
    while(b) {
        if( b & 1) res = (long long)res * a % mod;
        a = (long long)a * a % mod;
        b >>= 1;
    }
    return res;
}

int Inv(int x)
{
    return qpow(x, mod - 2);
}

int C(int a, int b)
{
    if(a < b) return 0;
    int res = 1;
    for(int i = 1, j = a; i <= b; ++i, -- j)
    {
        res = res * j % mod;
        res = res * inv[i] % mod;
    }
    return res;
}

int lucas(int a, int b, int p)
{
    if(a < p && b < p) return C(a, b);
    return C(a % p, b % p) * lucas(a / p, b / p, p);
}

signed main()
{
    for(int i = 1; i <= 20; ++i) inv[i] = Inv(i);
    // 预处理所有逆元（直接算超时）
    cin >> n >> m;
    int ans = 0;
    for(int i = 1; i <= n; ++i) cin >> a[i];
    for(int i = 0; i < 1 << n; ++i)
    {
        if(i == 0)
        {
            ans = (ans + lucas(n + m - 1, n - 1, mod)) % mod;
        }
        else 
        {
            long long t = n + m;
            int p = 0;
            for(int j = 0; j < n; ++j)
            {
                if(i >> j & 1 )
                {
                    t -= a[j + 1];
                    p ++;
                }
            }
            t -= p + 1;
            if(p & 1) ans = (ans - lucas(t, n - 1, mod)% mod + mod) % mod;
            else ans = (ans + lucas(t, n - 1, mod)) % mod;
        }
    }
    cout << (ans % mod + mod) % mod << endl;
}

Mobius函数

设正整数$N$按照算术基本定理分解质因数可以写成：$N=p_1^{c_1}{p}_2^{c_2}\cdots p_m^{c_m}$，定义函数
$$\mu (x)=\{\begin{array}{l}0\exists i\in [1,m]，c_i>1\\ 1m\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}2),\forall i\in [1,m],c_i=1\\ -1m\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}2),\forall i\in [1,m],c_i=1\end{array}$$
求某一项 $Mobius$的值，可以分解质因数计算.
若要求$1\sim n$的每一项$Mobius$的值。我们依旧可以利用埃式筛，或者线性筛来计算，时间复杂度接近$O(n)$
代码：

// 埃筛
void Mobius(int n)
{
    for(int i = 1; i <= n; ++i) mobius[i] = 1, st[i] = 0;
    for(int i = 2; i <= n; ++i)
    {
        if(!st[i])
        {
            mobius[i] = -1;
            for(int j = i + 1; j <= n; j += i)
            {
                st[j] = 1;
                if((j / i) % i == 0) mobius[j] = 0;
                mobius[j] *= -1;
            }
        }
    }
}
// 线性筛
void Mobius(int n)
{
	mobius[1] = 1;
    for(int i = 2; i <= n; ++ i)
    {
        if(!st[i])
        {
            prime[++ idx] = i;
            mobius[i] = -1;
        }
        for(int j = 1; i * prime[j] <= n; ++j)
        {
            st[i * prime[j]] = true;
            if(i % prime[j] == 0) {
                mobius[i * prime[j]] = 0;
                break;
            }
            mobius[i * prime[j]] = - mobius[i];
        }
    }
}

例题2：

题目链接

题意：

达达正在破解一段密码，他需要回答很多类似的问题：
对于给定的整数 $a,b$ 和 $d$，有多少正整数对 $x,y$，满足 $x\le a，y\le b$，并且 $gcd(x,y)=d$。
作为达达的同学，达达希望得到你的帮助。

Sol:

首先：$$gcd(a,b)=d\iff gcd(a/d,b/d)=1$$
证明很简单：$x,y$的最大公因数是$d$那么他们分别除了$d$以后一定是互质的。很显然
那么问题转化为：求$x\le a/d，y\le b/d$满足$gcd(x,y)=1$的个数。我们考虑容斥，
则答案为：所有的对数 - 不互质的对数
用式子表示：
$$\sum _{i=1}^{min(a,b)}\lfloor \frac{a}{i}\rfloor \ast \lfloor \frac{b}{i}\rfloor \ast \mu (i)$$
证明：
令：$a=\lfloor \frac{a}{d}\rfloor ,b=\lfloor \frac{b}{d}\rfloor$, $n=min(a,b)$，那么所有的对数：$a\ast b$

接下来要找出$gcd(x,y)>1$的数量。
不妨设集合$S_i,T_i$,分别表示$[1,a]$中因子有$i$的数的个数$\lfloor \frac{a}{i}\rfloor$，$[1,b]$中因子有$i$的数的个数$\lfloor \frac{b}{i}\rfloor$。那么容易知道，两个数不互质，那么他们的公因子大于等于$2$。
根据容斥定理和$Mobius$函数的定义，则有不互素的个数为:
$$\sum _{i=2}^n\lfloor \frac{a}{i}\rfloor \ast \lfloor \frac{b}{i}\rfloor \ast \mu (i)$$

此时，时间复杂度以及是$O(n)$，但此题仍然会超时。考虑分块做法，时间复杂度$O(\sqrt{n})$

Code：

#include <bits/stdc++.h>
using namespace std;
const int N = 50010;

int prime[N], mobius[N], sum[N], idx;
int n, a, b, d;
bool st[N];

void Mobius(int n)
{
    for(int i = 1; i <= n; ++i) mobius[i] = 1, st[i] = 0;
    for(int i = 2; i <= n; ++i)
    {
        if(!st[i])
        {
            mobius[i] = -1;
            for(int j = i + 1; j <= n; j += i)
            {
                st[j] = 1;
                if((j / i) % i == 0) mobius[j] = 0;
                mobius[j] *= -1;
            }
        }
    }
}

void init( int n )
{
    mobius[1] = 1;
    for(int i = 2; i <= n; ++ i)
    {
        if(!st[i])
        {
            prime[++ idx] = i;
            mobius[i] = -1;
        }
        for(int j = 1; i * prime[j] <= n; ++j)
        {
            st[i * prime[j]] = true;
            if(i % prime[j] == 0) {
                mobius[i * prime[j]] = 0;
                break;
            }
            mobius[i * prime[j]] = - mobius[i];
        }
    }
    for(int i = 1; i <= n; ++i) sum[i] = sum[i - 1] + mobius[i];
}

int main()
{
    init(N - 1);
    int tt;
    scanf("%d", &tt);
    while(tt -- )
    {
        scanf("%d%d%d", &a, &b, &d);
        a /= d, b /= d;
        n = min(a, b);
        long long ans = 0;
        for(int l = 1, r; l <= n; l = r + 1)
        {
            r = min(n, min(a / (a / l), b / (b / l)));
            ans += (sum[r] - sum[l - 1]) *(long long)(a / l) * (b / l);
        }
        printf("%lld\n", ans);
    }
}