抱歉，您的浏览器无法访问本站
本页面需要浏览器支持（启用）JavaScript
了解详情 >

数论补完计划 Part5 组合计数

oi算法

发布于：2024年9月25日

字数：2.9k 字

时长：11 分钟

排列数与组合数

排列组合基础

对于一个正整数 ，它的阶乘等于 ，记作 。特殊地，。

小学时我们就知道，用 三个数字最多能表示出 个互不相同的三位数字；若改为 组成四位数，那么答案就该是 种。这其实就是排列数的一个经典应用。探究这个答案是怎么得到的——对于第一个位置，可以有 种填数方法；对于第二个位置，因为先前已经用掉一个数了，此时可选的数就只剩 个，根据分步乘法原理，就应该乘上 ，以此类推……我们会惊奇地发现答案其实就是 ，或者可以表示成 。

排列数 是指从 个物体中选出 个并进行排列的总方案数。公式为 。排列所隐含的意思是，选中的 个物品之间是存在顺序差别的。例如上例选数字， 和 显然是不相同的两个数字，因此需要用到排列数。这也是它和组合数最大的差别。

上了高中，我们又知道，从 个人中选出 个人组队，总共能有 种组队方案。答案是 吗？鉴于 一队和 一对是等价的，那么其实队内的排序是不必要的，考虑到两个人能变出 种排序，意味着同一个方案会出现两次，答案其实是 。

组合数 （或者用二项式表示法表示为 ）是指从 个物体里选出 个并进行组合的总方案数。公式为 。它和排列数一样，只要 均为正整数，那么结果也都是整数。特殊地，。组合蕴含的意思是，不考虑选中元素之间的顺序关系，正如选择 组成一组和选择 成一组是完全等价的。

记不住组合数公式的看这张图：

还有特殊情况，例如 ，当 且 时，我们规定计算结果是 。

一个重要的恒等式是：。可以套公式验证。

排列组合扩展

比如说二项式定理，这个非常的好用，尤其是当你需要计算 的时候。二项式定理说道，对于 的展开形式，我们可以通过公式：

来快捷得到。注意 从 开始。 时我们很熟悉，完全平方展开，系数分别是 ； 还好，系数为 ； 时就是 。很容易发现这个系数似乎关于中间那一项左右对称，这一点可以用恒等式 完美解释。

跟二项式定理密切联系的还有一个杨辉三角，它长这个样子：

发现所有的二项式系数都可以在对应行找到。因此可以说，杨辉三角第 行的第 个数其实就对应组合数 。再次仔细观察，可以发现当前的数字可以由它肩上的两个数字相加得到，例如 ，于是所有的组合数都可以由这个递推式得到：。

杨辉三角还有许多有趣的性质，例如将第 行的所有二项式系数加起来能得到 ，要证明也很简单，利用二项式定理代入特值 ，这样就会得到 。还有，你还可以从第 行开始，将杨辉三角中的数按特定规则相加，就可以得到斐波那契数列的项，见下图：

即 。可以在了解组合数递推公式的基础上理解这个规律。

Lucas 定理

用于组合数取模问题，公式其实很简洁：

其中 为质数（暗示扩展），除法为整除。而当 时，第一项可以通过递归这个公式得到。

如果考虑递归，我们其实发现它和 进制分解蛮像。事实上， 定理就可以看作是 在 进制下进行对位求组合数的结果之积。也就是说，假如 的 进制表示分别为 和 （不存在对应为则为 ），那么 。证明略（不会）。它经常和数位DP放在一起考。

虽然它简化了计算，但是它本质上还是需要计算组合数的值。因此需要维护模意义下的阶乘和阶乘的逆元，或者你还可以暴力求解，适用于模数易变的情况。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// 暴力求解组合数
int C(int n, int m, int p) {
    int res = 1;
    for (int i = 1, j = n; i <= m; i++, j--) {
        res = ((ll) res % p * j % p * inv(i, p) % p) % p;
    }
    return res;
}

int lucas(int n, int m, int p) {
    if (m > n) return 0; // 必须特判，否则 RE+WA
    if (n < p && m < p) return C(n, m, p) % p;
    return (lucas(n / p, m / p, p) % p * lucas(n % p, m % p, p) % p) % p;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
// 预处理
int lucas(int n, int m, int p) {
    if (n < m) return 0; // 必须特判
    if (n < p && m < p) return (fac[n] % p * infac[m] % p * infac[n - m] % p) % p;
    return (lucas(n / p, m / p, p) % p * lucas(n % p, m % p, p) % p) % p;
}

void init(int lim) {
    fac[0] = infac[0] = 1;
    for (int i = 1; i <= lim; i++) {
        fac[i] = (fac[i - 1] * i % p) % p;
        infac[i] = (infac[i - 1] * inv(i, p) % p) % p;
    }
}

无论你选用哪一种，都不要忘了对特殊情况的特判，例如 的答案应为 。

扩展 Lucas 定理

回收伏笔

既然 定理只适用于 为质数，何不把模数 拆成质数呢？根据算术基本定理，这是可以做到的，假设最终分解为 。我们把代求的组合数当作一个未知数 ，联立可以得到方程组：

发现可以使用朴素的中国剩余定理合并（不同质数之间显然互质），前提是我们得知道 的值。

根据组合数的阶乘定义式得到：。但是还不能直接算，万一逆元不存在（ 或 ）怎么办？既然导致无解的原因是出现了 的倍数，那么我们就把 因子拆除了！设 分别为 中 因子的个数，那么原问题可以变为求解：

此时因为除去了因子 ，逆元就可以放心求解了。此时还面临一个问题，即形如 的算式的快速求解。分母可以用逆元求，其实时间复杂度瓶颈还是在阶乘取模上，因此只需要解决 即可。我们从一个实例来感受如何求解：

例 3.1

求解 的值

即 。我们先写阶乘，再把 的倍数提出来：

发现原式分为三个部分： 的 次幂、 的阶乘和与 互质的数的乘积。对于第二部分，由于可能仍然存在 的倍数，故递归求解；对于第三部分，有一个很好的性质——把互质的数按 个一组分下去（余下的另算），发现每个完整组的乘积模 的值是相同的，故快速幂解决，最后暴力乘上剩下的即可。

假设每一个完整组的乘积模 的值均为 ，剩余组模出来的值为 ，那么可以得到：

递归的下一层就是求解 的值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef pair<ll, ll> Equation;

ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    ll d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

ll inv(ll a, ll p) {
    ll ans, tmp;
    exgcd(a, p, ans, tmp);
    return (ans % p + p) % p;
}

ll qmul(ll a, ll b, ll p) {
    ll res = 0;
    while (b) {
        if (b & 1) res = (res + a) % p;
        a = (a + a) % p;
        b >>= 1;
    }
    return res % p;
}

ll qpow(ll a, ll b, ll p) {
    ll res = 1;
    while (b) {
        if (b & 1) res = (res % p * a % p) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return res % p;
}

ll mul(initializer_list<ll> args, ll p) {
    // 很方便的变长参数列表求连乘，int128 防止爆范围
    __int128 res = 1;
    for (ll i: args) {
        res = (res % p * i % p) % p;
    }
    return (ll) (res % p);
}

ll CRT(vector<Equation> &v) {
    // 中国剩余定理部分
    ll M = 1;
    ll res = 0;
    for (Equation eq: v) M *= eq.second;
    for (Equation i: v) {
        ll Mi = M / i.second;
        res = (res + qmul(qmul(i.first, Mi, M), inv(Mi, i.second), M)) % M;
    }
    return res % M;
}

ll resolve(ll a, ll p, ll ppow) {
    // 递归过程求解，即实现例3.1
    if (!a) return 1; // 1的阶乘模质数当然无论如何都是1
    ll s = 1;
    for (ll i = 1; i <= ppow; i++) {
        if (i % p) s = mul({s, i}, ppow); // 完整组循环节连乘
    }
    s = qpow(s, a / ppow, ppow); // 对完整组的积求幂
    for (ll i = a / ppow * ppow + 1; i <= a; i++) {
        if (i % p) s = mul({s, i}, ppow); // 对不完整组求积
    }
    return mul({s, resolve(a / p, p, ppow)}, ppow); // 递归过程
}

ll getA(ll n, ll m, ll p, ll ppow) {
    // 获取 CRT 方程组中的 a
    ll x = 0, y = 0, z = 0;
    // 统计 n, m, n - m 分别是 p 的多少次幂
    for (ll i = n; i; i /= p) x += i / p;
    for (ll i = m; i; i /= p) y += i / p;
    for (ll i = n - m; i; i /= p) z += i / p;
    // 套公式
    return mul({
        qpow(p, x - y - z, ppow),
        resolve(n, p, ppow),
        inv(resolve(m, p, ppow), ppow),
        inv(resolve(n - m, p, ppow), ppow)
    }, ppow);
}

ll exLucas(ll n, ll m, ll p) {
    // 分解质因数
    vector<Equation> E;
    for (ll i = 2; i * i <= p; i++) {
        if (p % i == 0) {
            // 统计对应质数以及它的幂
            ll c = 1;
            while (p % i == 0) {
                c *= i;
                p /= i;
            }
            E.emplace_back(getA(n, m, i, c), c);
        }
    }
    if (p > 1) E.emplace_back(getA(n, m, p, p), p); // 最后一个质数
    return CRT(E);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

    ll n, m, p;
    cin >> n >> m >> p;
    cout << exLucas(n, m, p) << endl;

    return 0;
}

其中 getA(n, m, p, ppow) 用来求解 ，即中国剩余定理中的 ； resolve(a, p, ppow) 实现求解 。

Lucas 定理的实际应用

鉴于 定理在进制分解上的性质，它其实可以跟一些涉及到进制的算法放在一起考。例如数位DP，数位DP相关内容请见：记忆化搜索 数位DP。

洛谷 P7976 [Stoi2033] 园游会

题目地址：P7976

题目难度：提高+/省选-

设 ，给定 ，求：

对 取模。

对于 的数据，。

根据 定理，我们知道 。那怎么把这个用到函数 上呢？我们可以发现， 满足 ，即积性函数。

我们把 转为三进制，那么此时求 就可以通过按位算组合数并连乘得到。固定一个 ，此时对于 的每一位，有几种可能：

  1. ，保证 的情况下，，此时贡献为
  2. ， 只能枚举 ，贡献为
  3. ，此时 ，贡献为

因此，真正能对答案产生实质性贡献的就是 的情况。换句话说，当 在 内枚举时，我们只需要统计 取不同值时（三进制下） 位的个数即可，也就是标准的数位DP求解 内数在三进制下总共有多少个 。记结果为 ，那么答案就是 。

洛谷 P6669/BZOJ 4737 [清华集训2016] 组合数问题

题目地址：P6669

题目难度：省选/NOI-

题目来源：清华集训  2016

给定 和 ，对于所有的 有多少对 满足 是 的倍数。答案对 取模。

对于 的测试点， ，，且 是一个质数。

读题，倍数必定满足关系 。根据 定理的基数对位组合求积意义，我们可以把 都转成 进制，并且，只要在连乘式中出现任意一个及以上的 项，整个结果就是 ，也即 的倍数。根据组合数的计算公式，发现只要 满足 ，那么结果就是 。

布尔型维护当前是否已经存在 项。此时考虑同时填两个不同数的相同位置，也就是一次性填 两个数，自然需要用到双层循环。鉴于题目要求 ，我们不仅需要维护 是否顶到上界 ，还要额外维护 是否顶到上界 和上界 。总的状态数是 ，转移是双层循环的 ，时间复杂度约为 ，完全足够。

题解同步于本站

本页文档最后更新于：2024年9月28日

算法

数论

oi

评论