容斥原理

容斥原理

1. 容斥原理

原理

• 容斥原理属于集合论的内容。
• 如果使用|S|表示集合S中元素的个数，那么高中我们最常见的是三个集合的情形：给定三个集合 $S_1、S_2、S_3$，求这三个集合元素的个数（重复的元素技能计算一次），如下图：

$$|S_1\cup S_2\cup S_3|=|S_1|+|S_2|+|S_3|-|S_1\cap S_2|-|S_1\cap S_3|-|S_2\cap S_3|+|S_1\cap S_2\cap S_3|$$

• 推广：对于n个集合，让求解一共出现了多少个元素，则有如下公式：

$$\left|\bigcup _i{S}_i\right|=\sum _i|S_i|-\sum _{i,j}|S_i\cap S_j|+\sum _{i,j,k}|S_i\cap S_j\cap S_k|-......$$

上式就是容斥原理内容。

• 如下是对容斥原理的证明：思路是查看每个元素是否是只被计算了一次，如果是的话，说明公式是正确的。

  不妨设元素x在k个集合（$1\le k\le n$）中出现过，根据公式则其被计算的次数为：
  $$C_k^1-C_k^2+C_k^3-...+(-1{)}^{k-1}{C}_k^k$$
  根据二项式定理，可知：
  $$C_k^1-C_k^2+C_k^3-...+(-1{)}^{k-1}{C}_k^k=-(1-1{)}^k+C_k^0=1$$
  因此每个元素只计算了一次，该公式是正确的。

• 另外我们考虑一个容斥原理中等式右侧有多少项相加减，第一项相当于从n个元素中拿1个，第二项相当于从n个元素中拿2个，…，因此相加减的项数有：
  $$C_n^1+C_n^2+C_n^3-...+C_n^n$$
  那么这个式子等于多少呢？如果上面再加上一项$C_n^0$，表示从n个元素中拿0个，则相当于问n个元素有多少个子集，考虑每个元素都可选可不选，因此一共有$2^n$个子集，因此有：
  $$C_n^0+C_n^1+C_n^2+C_n^3-...+C_n^n=2^n$$
  整理可以得到：
  $$C_n^1+C_n^2+C_n^3-...+C_n^n=2^n-1$$
  一共有$2^n-1$项相加减。

2. AcWing上的容斥原理题目

AcWing 890. 能被整除的数

问题描述

• 问题链接：AcWing 890. 能被整除的数

  

分析

• 最简单的思路是对于每个整数x，判断这m个质数中是否存在能整数x，如果存在，答案加一，这样做法的时间复杂度是$O(n\times m)$的，会TLE，不可取。

• 另一种思路是使用容斥原理，考察能被 $p_i(1\le i\le m)$，能被 $p_i$ 整除的数的数量为 $\lfloor \frac{n}{p_i}\rfloor$，结果加上这些数量，但是会有重复，比如同时能被两个质数 $p_i,p_j(1\le i,j\le m)$ 整除的数据被计算了两次，需要减去，…，因此最终的答案：

$$\sum _i\lfloor \frac{n}{p_i}\rfloor -\sum _{i,j}\lfloor \frac{n}{p_i\times p_j}\rfloor +\sum _{i,j,k}\lfloor \frac{n}{p_i\times p_j\times p_k}\rfloor -......$$

• 时间复杂度：一共有$2^m-1$项相加减，每次计算主要是质数相乘，最多相乘m-1次，因此时间复杂度为$O(2^m\times m)$。

• 另外还要考虑如何遍历所有的可能组合，这有一个常用的技巧，我们可以通过二进制的方式进行遍历，即从1循环到$2^m-1$，对于其中的每个数据，对应一种集合，该数据的二进制表示中如果某一位是1，代表需要乘上该质数。

代码

• C++

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 20;

int p[N];  // 输入的质数

int main() {
    
    int n, m;
    cin >> n >> m;
    for (int i = 0; i < m; i++) cin >> p[i];
    
    int res = 0;
    for (int i = 1; i < 1 << m; i++) {
        
        int t = 1, s = 0;  // t表示当前i对应集合中所有质数乘积, s表示集合中数据个数
        for (int j = 0; j < m; j++) 
            if (i >> j & 1) {
                if ((LL)t * p[j] > n) {
                    t = -1;  // 乘积大于n，res不用变
                    break;
                }
                t *= p[j];
                s++;
            }
        
        if (t != -1) {
            if (s % 2) res += n / t;
            else res -= n / t;
        }
    }
    
    cout << res << endl;
    
    return 0;
}

AcWing 214. Devu和鲜花

问题描述

• 问题链接：AcWing 214. Devu和鲜花

  

分析

• 首先我们考虑一个简化版的该问题，假设每个盒子中的花有无限多个，我们有多少中种选择方案？

• 假设每个盒子中选择的花的数量为$x_i$，则$x_1+x_2+...+x_N=M,其中x_i\ge 0$，我们令$y_i=x_i+1$，则$y_1+y_2+...+y_N=M+N,其中y_i>0$，我们可以使用隔板法解决该问题，类似于组合计数中的AcWing 1312. 序列统计，不过本题是等号，AcWing1312是小于等于。

• M+N个小球可以在M+N-1个缝隙中放置N-1个隔板，因此存在$C_{M+N-1}^{N-1}$种方案。

• 然后我们考虑限制条件，这里要求$x_i\le A_i$，利用补集的思想，我们考虑不满足各个限制条件的情况，让$S_i$表示不满足$x_i\le A_i$的所有方案的集合（即在这个集合中一定有$x_i>A_i$），例如$S_1$表示不满足第一个条件的集合（其他条件满不满足无所谓）。则总体的方案数（所有方案-不合法的方案）：

$$C_{M+N-1}^{N-1}-|S_1\cup S_2\cup ...\cup S_N|$$

• 其中不合法的方案数可以使用容斥原理求解。考虑$|S_1\cup S_2\cup ...\cup S_N|$的求解，将该式展开得到：

$$\left|\bigcup _i{S}_i\right|=\sum _i|S_i|-\sum _{i<j}|S_i\cap S_j|+\sum _{i<j<k}|S_i\cap S_j\cap S_k|-......$$

• 我们考虑$|S_1|$的求解，根据定义，$S_1$对应的集合中一定有$x_1>A_1$，我们可以首先让$x_1=A_1+1$，之后可以在N个盒子中选取$M-(A_1+1)$只花（每个盒子中的花可以看成无限的），这就是$S_1$对应的所有方案；从N个盒子中选取$M-(A_1+1)$只花（每个盒子中的花可以看成无限的）就是我们刚开始考虑的问题，因此$S_1$对应的方案数为$C_{M-(A_1+1)+N-1}^{N-1}$。

• 类似的考虑$|S_1\cap S_2|$的计算，则我们需要让$x_1=A_1+1，x_2=A_2+1$，之后可以在N个盒子中选取$M-(A_1+1)-(A_2+1)$只花（每个盒子中的花可以看成无限的），这就是$|S_1\cap S_2|$对应的所有方案；因此$|S_1\cap S_2|$对应的方案数为$C_{M-(A_1+1)-(A_2+1)+N-1}^{N-1}$。

• 以此类推我们可以求出不合法的方案数量。因此我们最终的合法方案数量为：

$$res=C_{M+N-1}^{N-1}-\sum _i{C}_{M+N-A_i-2}^{N-1}+\sum _{i<j}{C}_{M+N-A_i-A_j-3}^{N-1}-......$$

• 考虑上式如何求解，容斥原理展开式中一共有$2^N-1$项，上式中第一项可以看成一个$S$也不选的情况，一共有$2^N$项，可以使用二进制数据枚举所有情况，因此本题的时间复杂度为$2^N\times N$，这里乘以N是因为我们求解$C_{M+N-1}^{N-1}$的计算量为$O(N)$，因为N很小，我们按照定理求解$C_{M+N-1}^{N-1}$即可，即：

$$C_{M+N-1}^{N-1}=\frac{(M+N-1)\times (M+N-2)\times ...\times (M+1)}{(N-1)!}$$

• 我们使用up、down分别表示上式的分子和分母的逆元（可以使用快速幂求解，因为两者互质且题目给的模数为质数），我们发现分母是不会变的，因此可以只计算一次，这样本题就不会存在最后一个数据无法通过的情况了。

代码

• C++

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 25, mod = 1e9 + 7;

LL n, m;
LL A[N];
int down = 1;

int qmi(int a, int k, int p) {
    int res = 1 % p;
    while (k) {
        if (k & 1) res = (LL)res * a % p;
        a = (LL)a * a % p;
        k >>= 1;
    }
    return res;
}

int C(LL a, LL b) {
    if (a < b) return 0;
    int up = 1;
    for (LL i = a; i > a - b; i--) up = i % mod * up % mod;  // i本身就是LL
    
    return (LL)up * down % mod;
}

int main() {
    
    cin >> n >> m;
    for (int i = 0; i < n; i++) cin >> A[i];
    
    // 计算(n-1)!对mod的逆元
    for (int i = 1; i <= n - 1; i++) down = (LL)down * i % mod;
    down = qmi(down, mod - 2, mod);
    
    int res = 0;
    for (int i = 0; i < 1 << n; i++) {
        LL a = m + n - 1, b = n - 1;
        int sign = 1;  // 选择了偶数项S系数为正，否则为负
        for (int j = 0; j < n; j++) 
            if (i >> j & 1) {
                sign *= -1;
                a -= A[j] + 1;
            }
        res = (res + C(a, b) * sign) % mod;
    }
    
    cout << (res + mod) % mod << endl;
    
    return 0;
}

AcWing 215. 破译密码

问题描述

• 问题链接：AcWing 215. 破译密码

  

分析

• 本题中的需要用到知识点：莫比乌斯函数（Mobius函数）。下面给出该函数的定义：

• 对于给定正整数x，如果x可以被质因数分解为$p_1^{{\alpha }_1}\ast p_2^{{\alpha }_2}\ast ...\ast p_k^{{\alpha }_k}$，莫比乌斯函数记作$\mu (x)$，则：

$$\mu (x)=\{\begin{array}{l}0至少存在一个\alpha 的值大于1\\ 1所有\alpha 值都为1，且k为偶数\\ -1所有\alpha 值都为1，且k为奇数\end{array}$$

• 那么该函数有什么作用以及如何求解呢？

• 作用如下（存疑？）：可以求解1~N中和N互质的数的个数，我们使用$S_{p,q}$表示表示1~N中和N存在公因数p和q的所有数的集合（其中p、q都是质数），一共有$\lfloor \frac{N}{p}\rfloor$个数。则1~N中和N互质的数的个数为：

$$N-|S_2|-|S_3|-|S_5|-...+|S_{2,3}|+|S_{3,5}|+...-|S_{2,3,5}|-...$$

观察可以发现，上式中除了N之外的的其他项的系数就是莫比乌斯函数，因此上式可以改写成：
$$N-\sum \mu (p)\times |S_p|$$

• 如何求解1~N中每个数对应的莫比乌斯函数值呢？可以在质数线性筛法中求解，代码以及分析如下：

// 线性选法，时间复杂度：O(n)
void get_primes(int n) {
    mobius[1] = 0;
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes[cnt++] = i;
            mobius[i] = -1;  // i本身就是质数，k为奇数
        }
        for (int j = 0; primes[j] <= n / i; j++) {
            int t = primes[j] * i;
            st[t] = true;
            if (i % primes[j] == 0) {
                mobius[t] = 0;  // t至少存在两个质因子primes[j]
                break;
            }
            // mobius[i] = 0 时, mobius[t]也为0
            // mobius[i] != 0 时, 此时t比i多一个i质因数分解中不存在的质数, mobius[t]和mobius[i]异号
            mobius[t] = mobius[i] * -1;
        }
    }
}

• 现在我们考虑本题如何求解，问题是给定我们a、b、d，我们需要计算数对的数量，这些数对(x, y)需要满足：$1\le x\le a,1\le y\le b$且gcd(x, y)==d。

• 上述问题等价于：$1\le x^{\prime }\le \lfloor \frac{a}{d}\rfloor ,1\le y^{\prime }\le \lfloor \frac{b}{d}\rfloor$且gcd(x', y')==1的数对的数量。证明也很简单，只需要做一个简单的映射即可，即$x^{\prime }=\frac{x}{d},y^{\prime }=\frac{y}{d}$。

• 首先让$a=\lfloor \frac{a}{d}\rfloor ,b=\lfloor \frac{b}{d}\rfloor$。

• 现在问题转成了求互质的数的对数，这两个数的范围分别是从1~a，和从1~b。可以使用容斥原理。最后的答案是：

$$a\times b-\lfloor \frac{a}{2}\rfloor \times \lfloor \frac{b}{2}\rfloor -\lfloor \frac{a}{3}\rfloor \times \lfloor \frac{b}{3}\rfloor +\lfloor \frac{a}{6}\rfloor \times \lfloor \frac{b}{6}\rfloor +......$$

其中$\lfloor \frac{a}{2}\rfloor \times \lfloor \frac{b}{2}\rfloor$代表a、b存在公因数2的数的对数。上式可以形式化的写成：
$$a\times b+\sum _{i=1}^{min(a,b)}\mu (i)\times \lfloor \frac{a}{i}\rfloor \times \lfloor \frac{b}{i}\rfloor$$

• 根据上式我们可以看出每次查询的时间复杂度是$O(n)$的，因此总体的时间复杂度为$O(n^2)$，n最大为50000，会超时。因此需要优化。

• 我们考虑对于给定正整数n，考虑$\lfloor \frac{n}{1}\rfloor ,\lfloor \frac{n}{2}\rfloor ,\lfloor \frac{n}{3}\rfloor ,...,\lfloor \frac{n}{n}\rfloor$中实际存在多少不同的数据，另外如何求解出这些数据的个数。

• 存在多少不同的数据：

• 实际上上述数据的个数不超过$2\times \sqrt{n}$个，这是因为我们可以将数据分为两组：

  （1）$\lfloor \frac{n}{1}\rfloor ,\lfloor \frac{n}{2}\rfloor ,...,\lfloor \frac{n}{\sqrt{n}}\rfloor$；

  （2）$\lfloor \frac{n}{\sqrt{n}+1}\rfloor ,\lfloor \frac{n}{\sqrt{n}+2}\rfloor ,...,\lfloor \frac{n}{n}\rfloor$。

• 第（1）部分只有$\sqrt{n}$项；第（2）部分虽然项数很多，但是数值都是小于$\sqrt{n}$的，因此最多有$\sqrt{n}$项；因此不同的数据的个数不超过$2\times \sqrt{n}$个。

• 如何求解出这些数据的个数：

• 我们首先考虑如下问题：对于n，给定我们一个小于等于n的正整数x，我们希望找到值等于$\lfloor \frac{n}{x}\rfloor$且最大的一个数y满足$\lfloor \frac{n}{x}\rfloor =\lfloor \frac{n}{y}\rfloor$，我们可已看出y是x的函数，令$g(x)=y$，例如给定n=24, x=9，则g(x)=12，因为$\lfloor \frac{24}{9}\rfloor =\lfloor \frac{24}{12}\rfloor$，且12是最大的一个。

• g(x)的形式化定义如下：$\lfloor \frac{n}{x}\rfloor =\lfloor \frac{n}{g(x)}\rfloor$且$\lfloor \frac{n}{x}\rfloor >\lfloor \frac{n}{g(x)+1}\rfloor$。这里的g(x)是由公式解的，如下：

$$g(x)=\lfloor \frac{n}{\lfloor \frac{n}{x}\rfloor }\rfloor$$

• 我们可以证明g(x)是满足形式化定义中的两个式子的，首先证明：$\lfloor \frac{n}{x}\rfloor =\lfloor \frac{n}{g(x)}\rfloor$。（证明方式：A>=B, A<=B，则A==B）
• 根据g(x)的表达式，我们可以推出：

$$g(x)=\lfloor \frac{n}{\lfloor \frac{n}{x}\rfloor }\rfloor >=\lfloor \frac{n}{\frac{n}{x}}\rfloor =x$$

因此有：
$$\lfloor \frac{n}{x}\rfloor <=\lfloor \frac{n}{g(x)}\rfloor (1)$$
将g(x)带入$\lfloor \frac{n}{g(x)}\rfloor$，可以得到：
$$\lfloor \frac{n}{g(x)}\rfloor =\lfloor \frac{n}{\lfloor \frac{n}{\lfloor \frac{n}{x}\rfloor }\rfloor }\rfloor >=\lfloor \frac{n}{\frac{n}{\lfloor \frac{n}{x}\rfloor }}\rfloor =\lfloor \frac{n}{x}\rfloor (2)$$
因此（1）（2）式可知：$\lfloor \frac{n}{x}\rfloor =\lfloor \frac{n}{g(x)}\rfloor$。

• 下面接着证明形式化定义的第二个式子：$\lfloor \frac{n}{x}\rfloor >\lfloor \frac{n}{g(x)+1}\rfloor$。
• 令$n=k\times x+r,0\le r<x$（对x求带余除法），则等价于证明下式（此时$g(x)=\lfloor \frac{n}{k}\rfloor$）：

$$k>\lfloor \frac{n}{\lfloor \frac{n}{k}\rfloor +1}\rfloor$$

等价于证明：
$$k\times (\lfloor \frac{n}{k}\rfloor +1)>n$$
令$n=p\times k+q,0\le q<k$（对k求带余除法），带入上式可以得到：
$$k\times (p+1)>p\times k+q⟺k>q$$
因此第二个式子也是成立的。

---

• 有了上述理论之后，我们就可以快速求解出刚开始推出来的式子：

$$a\times b+\sum _{i=1}^{min(a,b)}\mu (i)\times \lfloor \frac{a}{i}\rfloor \times \lfloor \frac{b}{i}\rfloor$$

• $\lfloor \frac{a}{i}\rfloor \times \lfloor \frac{b}{i}\rfloor$对应的序列有很多值都是相同的，计算一次就行，假设值记为t，对应的区间假设为[l, r]，然后让t乘以$\mu (l)+\mu (l+1)+...+\mu (r)$即可求出该段中的和，$\mu (l)+\mu (l+1)+...+\mu (r)$可以使用前缀和求解。

• 时间复杂度为$O(n\times \sqrt{n})$，n次询问，每次询问有不多于$2\times \sqrt{n}$段数据相加。

代码

• C++

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 50010;

int primes[N], cnt;
bool st[N];
int mobius[N], sum[N];  // sum是mobius的前缀和

// 线性筛法求莫比乌斯函数
void init(int n) {
    mobius[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes[cnt++] = i;
            mobius[i] = -1;
        }
        for (int j = 0; primes[j] * i <= n; j++) {
            int t = primes[j] * i;
            st[t] = true;
            if (i % primes[j] == 0) {
                mobius[t] = 0;
                break;
            }
            mobius[t] = mobius[i] * -1;
        }
    }
    
    for (int i = 1; i <= n; i++) sum[i] = sum[i - 1] + mobius[i];
}

int main() {
    
    init(N - 1);
    
    int T;
    scanf("%d", &T);
    while (T--) {
        int a, b, d;
        scanf("%d%d%d", &a, &b, &d);
        a /= d, b /= d;
        int n = min(a, b);
        LL res = 0;
        for (int l = 1, r; l <= n; l = r + 1) {
            r = min(n, min(a / (a / l), b / (b / l)));  // 下次跳到的位置应该是r+1
            res += (sum[r] - sum[l - 1]) * (LL) (a / l) * (b / l);
        }
        printf("%lld\n", res);
    }
    
    return 0;
}