组合数求模

0. 写在前面

在程序设计中，可能会碰到多种类型的计数问题，其中不少涉及到组合数的计算，所以笔者写下这么一篇文章，期望能解决一些常规的组合数求模问题。以下部分内容改编自AekdyCoin的《组合数求模》，而且为了感谢他对（懵懂的）笔者的启发，这篇文章的标题与其文章相同。另外，感谢Picks将多项式运算的技巧在中国进行推广，感谢51nod提供了许多有趣的数论题目，感谢fotile96开源了他的FFT模板。

在接下来的内容中，文章将围绕 C(n,m)modp 进行分析，其中 C(n,m) 表示从 n 个人中选出 m 个人跪sd0061的方案数。

author: skywalkert 
original article: http://blog.csdn.net/skywalkert/article/details/52553048 
last update time : 2017-03-30

1. 枚举法

这一部分的内容主要是利用基本定义

C(n,m)=n!m!(n−m)!=n(n−1)⋯(n−m+1)1⋅2⋯m

来进行分析的。

1.1 直接枚举

如果 n−m<m ，不妨将 n−m 视作 m ，所以我们只需要考虑 m≤n2 的情况。

若 1,2,⋯,m 在模 p 意义下有乘法逆元，则考虑直接计算分子、分母在模意义下的值，利用扩展欧几里得算法 O(logp) 求逆即可，这样做求 C(n,m) 的复杂度是 O(m) 的。注意 p 可以不是质数。

如果需要求 C(n,1),C(n,2),⋯,C(n,m) ，也可以利用上述方法做到 O(mlogp) ，这是因为

C(n,m)=C(n,m−1)⋅n−m+1m

如果能提前 O(m) 预处理 1,2,⋯,m 的逆元，则上述问题可以做到 O(m) ，这是因为

m=x⋅⌊mx⌋+(mmodx)⇒x⋅⌊mx⌋+(mmodx)≡0(modm)⇒x−1≡−⌊mx⌋(mmodx)−1(modm)

注意，当 m 不为质数时， (mmodx) 不一定与 m 互质，所以上述做法需要有 m 是质数的条件。

1.2 预处理

由于

C(n,m)=n(n−1)⋯(n−m+1)1⋅2⋯m=(n−1)⋯(n−m+1)1⋅2⋯(m−1)(1+n−mm)=C(n−1,m−1)+C(n−1,m),

所以能够 O(nm) 计算出 C(n,m) 。

注意到上述式子只涉及加法，所以不需要考虑是否有逆元的问题。

这种做法可以预先处理出一定范围内的所有组合数，将计算多个组合数转化为查表，适用于多次查询组合数的问题。

若 1,2,⋯,n 在模 p 意义下有乘法逆元，也可以预处理 n! 和 (n!)−1 ，当逆元可以线性预处理时（ p 为质数），阶乘和阶乘的逆元也可以 O(n) 预处理。

注意，当 n 远大于 m 时，或许直接枚举的做法会更好。

1.3 质因数分解

若 1,2,⋯,m 中存在数字在模 p 意义下没有逆元，则需要避免使用相应的除法，而是在运算前将分子和分母中的相关项相消，一个显而易见的想法便是将分子和分母所对应的阶乘化为质因数分解式，在质因数的幂次上进行加减法。计算一个数的幂可以利用快速模幂算法做到幂次的对数复杂度。

如果不考虑模 p 意义，则只需要考虑分母中所含的质因数，分子也只需要对这些质因数进行分解，其他剩余的部分可以当作整体，这样做是 O(mlogm) 的，可在分母超出存储范围但是精确值不超过存储范围时使用。

而考虑模 p 意义，则注意到需要消除的数都是与 p 的最大公约数大于1的，因此只需要考虑 p 的质因数，将与 p 互质的部分直接计算（互质的分母部分直接求逆元），不互质的部分维护成质因数的幂次即可，这样做是 O(mlogp) 的，与1.1同理可以求 C(n,1),C(n,2),⋯,C(n,m) 。

当 n 与 m 同阶时，可以利用线性的筛法维护一些额外信息，从而 O(n) 完成组合数的质因数分解。具体做法如下：

1. 用筛法计算出正整数 x (x≤n) 的一个质因数 prx 。
2. 维护一个数组 c1,c2,⋯,cn 使得 C(n,m)=∏ni=1ici ，初始化每个位置为0，然后使 c1,c2,⋯,cn 各加一， c1,c2,⋯,cm 各减一， c1,c2,⋯,cn−m 各减一。
3. 降序枚举 x (2≤x≤n) ，如果 x 不为质数，则将 xcx 拆分为 prxcx⋅(xprx)cx ，具体操作是修改 cx,cprx,cxprx 。
4. 此时 cx 不为0则必然 x 为质数，而且 x 为质数时 cx 非负，即为组合数的质因数分解。

当然，如果只有少量数字没有逆元，还可以考虑扩大模数。例如，只有 x 没有逆元，但是要除以一个 x ，则可以计算原式在模 x⋅p 意义下的值，再除以 x ，而如果有多个数字时，也可将模数扩大它们的乘积倍。但一般情况下，这种方法的效率不高。

2. 定理法

这一部分的内容将涉及初等数论中两个常用的定理，利用定理将复杂的问题简单化，从而便于分析。

2.1 卢卡斯定理 (Lucas’ theorem)

在数论中，卢卡斯定理将组合数 C(n,m) 在模质数 p 意义下的值表示成一系列组合数 C(ni,mi) 的乘积，其中 ni,mi 分别表示 n,m 在 p 进制下的第 i 位数字，换句话说，令 n=∑ki=0ni⋅pi,m=∑ki=0mi⋅pi ，则有

C(n,m)≡∏i=0kC(ni,mi)(modp)

其中，如果 ni<mi 则 C(ni,mi)=0 。另外，这个定理的证明还用到了一个常用于化简式子的公式：当 p 是质数时， (A+B)p≡Ap+Bp(modp) 。

当 n 远大于 p 时，利用卢卡斯定理，我们可以将组合数化简到小于 p 的情况，从而在套用1.2中的预处理阶乘和阶乘逆元时，复杂度从 O(n) 降到了 O(p) 。

假设能够轻松地（进制转换）得到 p 进制下的 n,m ，则只需要计算 O(k)=O(logn) 个组合数，所以复杂度是 O(p+logn) 。但是希望大家能注意到卢卡斯定理的局限性：它只适用于模质数的情况。

2.2 中国剩余定理 (Chinese Remainder Theorem)

中国剩余定理适用于化简模线性方程组，具体来说，如果有

⎧⎩⎨⎪⎪⎪⎪⎪⎪x≡r1(modm1)x≡r2(modm2)⋯x≡rk(modmk)

，其中 gcd(mi,mj)=1 (i≠j) ，那么利用中国剩余定理可以得到

x≡∑i=1kriM′iMi(modM)

与原方程组等价，其中 M=∏ki=1mi,Mi=Mmi ，而 M′i≡Mi−1(modmi) 。

中国剩余定理实际上利用 k 次求逆元将方程组进行了合并，而一般情况下可能出现方程组模数不互质的情况，有两种解决办法。一种是将模数分解质因数，划分成多个小的方程组，同质数的方程组再根据幂次判断合法性并合并。另一种是利用二元一次不定方程的通解判断合法性和合并方程组。一般来说，前者需要增加代码量，后者需要一定的推导。

利用中国剩余定理，我们可以将模 p 的问题转化为模其质因子幂次的问题。当 p 是square-free number（因子中的平方数只有 1）时，可以结合卢卡斯定理分别解决问题，再将结果合并。当 p 不是square-free number时，需要使用下文的一些方法。

3. 其他化归法

这一部分的内容将涉及笔者认为常用的其他化简、归纳方法。笔者不才，私以为将这些方法与前面的结论相结合，可以解决目前大部分组合数求模的问题。

3.1 阶乘表示法

与1.3类似，在这一节我们考虑将组合数表示成 p 的质因数的幂次和其他部分，尝试简化问题，但不同的是，由于中国剩余定理的支持，这一节我们只考虑 p 是质因数幂次的情况，即 p=qe 。

组合数可以表示成阶乘和阶乘逆元的乘积，所以我们只要考虑如何表示阶乘。设 n!=qkn⋅vn ，其中 vn 与 q 互质，则需要计算 kn 的精确值与 vnmodqe 。

阶乘 n! 可以表示成前 n 个连续正整数 1,2,⋯,n 的乘积，将这 n 个数中与 q 互质的数提取出来，这一部分的数乘积累积给 vn 。而不互质的数必然是 q,2⋅q,⋯,⌊nq⌋⋅q ，每个数字提取一个 q ，则至少可以提取出 ⌊nq⌋ 个 q 累加给 kn 。那么剩下没考虑的数字就是 1,2,⋯,⌊nq⌋ ，也即 ⌊nq⌋! ，再次重复上述过程进行提取即可。

提取 q 的过程是十分简单的，只需要 O(logn) 次提取。而主要的复杂度在于计算 1,2,⋯,n 中与 q 互质的数字乘积在模 qe 意义下的取值。由于 1,2,⋯,qe 和 qe+1,qe+2,⋯,2⋅qe 在模意义下对应相等，因此可以将 1,2,⋯,n 表示成 ⌊nqe⌋ 个 1,2,⋯,qe 和 1,2,⋯,nmodqe ，因此可以预处理前 x (x<qe) 个正整数里与 q 互质的数字乘积，再利用快速模幂算法高斯泛化的威尔逊定理处理前者，这样做的整体复杂度是 O(qe+logn) 。这个做法在 e=1 时等同于2.1中的预处理。

3.2 分治与多项式

继续分析上一节中的互质的数字乘积问题，在上一节中我们只用到了模 qe 的循环节，但是在这样的循环节中还是有 qe−1 个位置是不产生贡献的，剩下产生贡献的位置被分成了许多个长度为 q−1 的区间，当 e>1 时是十分浪费的。

考虑多项式 f(x)=(x+1)(x+2)⋯(x+q−1) ， f(0) 就表示了区间 [1,q) 中正整数的乘积， f(q) 就表示了 [q+1,2⋅q) 的情况，以此类推我们可以表示出循环节中每个长度为 q−1 的区间乘积，而前 n (n<qe) 个正整数里与 q 互质的数字乘积就可以表示为多项式在前 ⌊nq⌋ 个 q 的非负倍数处的取值乘积，再乘上不超过 q 个零碎的数字，问题再一次化简。

注意到 x 的取值是 q 的倍数时，这个多项式里 xi 的项的取值一定是 qi 的倍数，故当 i≥e 时 xi 的项在模意义下一定为 0 ，因此多项式可以在模 xe 意义下进行运算。（更高效而且更优美的性质是，多项式里  xi  的系数只需要在模  qe−i  意义下维护即可）

考虑这样一个事实：如果已知 [x+1,x+aq) 对应的多项式为 g(x) ，则 [x+1,x+2aq) 对应的多项式为 g(x)⋅g(x+aq) 。同理可以用 [x+1,x+aq) 和 [x+1,x+bq) 的结果得到 [x+1,x+(a+b)q) 的结果。而我们的运算都是在前 e 项进行的，因此利用类似快速幂的方法可以得到 [x+1,x+⌊nq⌋q) 对应的多项式，再求其常数项（即 f(0) ）与零碎的数字乘积的值即可，这样做是 O((q+e3)logn) 的，如果将  O(e)  次计算的顺序调整一下，复杂度可以做到  O((q+e2)logn)  ，并且不会使用到高斯泛化的威尔逊定理。

当 q≥e2>1 时， q=p1e=O(p√) ，总复杂度也可以写为 O(p√logn) ，可以完成一些相对来说比较大的情况。

3.3 分块与多项式

还是继续分析上一节中的互质的数字乘积问题，在上一节中我们对 e>1 的情况取得了较好的做法，在这一节中我们会对 e=1 的情况进行另外一个方向的分析，即组合数模大质数问题。

我们需要考虑的是计算 n! (n<p) 在模质数 p 意义下的值，直接求值会很慢，但直接构建多项式会使得构建过慢、求值很快。为了均衡二者，可以设定一个阈值 T ，构建多项式 Q(x)=∏Ti=1(x+i) ，尝试求出 Q(0),Q(1),⋯,Q(⌊n−T+1T⌋) ，与剩下不超过 T 个零碎的数字直接乘起来，从而得到 n! 。

多项式 Q(x) 是一个次数为 T 的多项式，利用分治+FFT构建这个多项式的复杂度是 O(Tlog2T) 。由于多项式 Q(x) 在一个点 x=xi 的取值即 Q(x)mod(x−xi) ，因此在多点求值可以利用分治的技巧，将求值的序列一分为二，同时将多项式降次成两个次数更小的多项式，这其中需要预处理构建一些辅助多项式，从而做到复杂度为 O(nTlog2nT) 。综上所述， T 的最优值应取 O(n√) ，限于常数问题可能需要自行调参，总体的复杂度为 O(n√log2n) 。

在上述过程中涉及到的系数运算都是模 p 意义下的，但直接使用FFT会使得过程中数值达到 O(Tp2) 的级别，利用修正许多误差且使用优化技巧的FFT可以解决这个问题。

3.4 更复杂的情况

如果对扩展Lucas定理感兴趣，可以参考Andrew Granville在1997年发表的文章Arithmetic Properties of Binomial Coefficients I: Binomial coefficients modulo prime powers，适合  e>1  的情况。但是笔者从论文中了解到的做法比3.2的做法要差一点，做法类似于拉格朗日插值法求出互质的数字乘积，而3.2中维护多项式的做法会更加灵活、高效。

可以发现，当 e>1 时3.2的做法比3.3的做法更好，但也不难发现，3.3中提到的分块技巧是可以与3.2中的分治技巧进行结合的（弥补其中  O(q) 的部分），如果对于更加一般化的做法感兴趣，欢迎与我联系。

4. 例题讲解

这一部分的内容将按照上文内容的顺序依次讲解适用相关方法的例题，目的是加深对相关方法的理解与应用。

4.1 Ural 1114 Boxes

简要题意

有 n 个盒子、 a 个红球和 b 个蓝球，每个盒子可以放任意多个红球和任意多个蓝球，盒子可以空着，球也可以不全放进去，问有多少种放求的方案数。

1≤n≤20,0≤a,b≤15

简要题解

红蓝球互不影响，分别计算方案数后相乘即可。而 m 个球放到 n 个盒子的方案数即 C(m+n−1,n−1) ，假设多一个盒子放置那些没放进去的球，则可知方案数为 C(m+n,n) ，即答案为 C(a+n,n)⋅C(b+n,n) 。初步估算最大值为 C(35,20)2≈1019 ，因此需要使用64位无符号整型数存储，为了不使运算过程中出现溢出，需要使用质因数分解的方法，套用1.2使用加法预处理组合数即可。

4.2 UVaLive 7040 color

简要题意

有 n 个连续的格子和 m 种不同颜色的染料，现在要用恰好 k 种染料对这 n 个格子染色，每种染料至少染一个格子，并且任意相邻的格子颜色不同，问这样选出染料并染色的合法方案数对 109+7 取模的值。

1≤n,m≤109,1≤k≤min(106,n,m)

简要题解

选出 k 种颜色的方案是 C(m,k) ，最后乘上即可。如果只是使得任意相邻的格子颜色不同，那么方案数是 k⋅(k−1)n−1 。现在要使得每种染料至少染一个格子，则可以考虑上述方案中只出现 i(i≤k) 个颜色的方案 f(i) ，则有 k⋅(k−1)n−1=∑i≤kC(k,i)⋅f(i) ，利用二项式反演可以得到 f(k)=∑i≤k(−1)k−iC(k,i)⋅i⋅(i−1)n−1 。因此只需要预处理逆元后直接枚举 C(k,i) 即可，时间复杂度 O(k+klnklogn) 。

4.3 COGS 1284 [HNOI2004] 树的计数

简要题意

已知一个有 n 个结点的树第 i 个结点的度数为 di ，问满足这样的条件的不同的树有多少棵。

n≤150 ，保证答案不超过 1018

简要题解

由树的Prüfer序列可知，第 i 个结点会在序列中出现 di−1 次，总共有 n−2 个数字，所以答案为 (n−2)!∏ni=1(di−1)! ，由于保证了答案不超过 1018 ，保险起见最好将式子表示成一些数字的乘积，因此需要套用1.3中的质因数分解方法。

4.4 BNUOJ 51637 Quailty’s Life

简要题意

从 (0,0) 出发走到 x=n 或 y=m ，每一步需要使横坐标或纵坐标加一，也可使二者同时加一，路上不能走到 k 个定点 (ai,bi) (i=1,2,⋯,k) ，问路径的方案数对 p 取模的值。

1≤n≤109,1≤m≤104,0≤k≤10,1≤p<231

简要题解

详细题解见BNUOJ 51637 Quailty’s Life | BZOJ 4587 推箱子。

首先对于 k 个定点进行去重，不难发现一个定点能到达另一个定点这样的关系构成一个有向无环图，因此可以尝试进行容斥，设 f(i) 表示从起点出发走到第 i 个定点且之前没有走到其他定点的方案数， way(x,y) 表示横坐标增加 x 、纵坐标增加 y 的路径方案数，则有 f(i)=way(ai,bi)−∑j≠if(j)⋅way(ai−aj,bi−bj) ，这只需要 O(k2) 次计算 way(x,y) 。

问题在于计算 way(i,j) ，而且对于最后的终止情况还需要计算 i=n 和 j=m 的情况之和。设从 (0,0) 走到 (i,j) 的一种方案使用了 a 次 x+1 ， b 次 y+1 ， c 次 x+1 同时 y+1 ，则有 

⎧⎩⎨a+c=ib+c=ja,b,c≥0⇒⎧⎩⎨a=i−cb=j−c0≤c≤min(i,j)

，对应的方案数是

(i+j−c)!(i−c)!(j−c)!c!=C(i+j−c,i)⋅C(i,c)=C(j+i−c,j)⋅C(j,c)

。对于最后的求和部分，先考虑 i=n 的情况，则 0≤j≤m 的情况为

∑j=0m∑c=0min(n,j)C(n+j−c,n)⋅C(n,c)=∑c=0min(n,m)C(n,c)∑j=cmC(n+j−c,n)=∑c=0min(n,m)C(n,c)⋅C(n+1−c,n+1),

那么也是枚举 c 就可以了，不用枚举 j 。对于 0≤i≤n,j=m 的情况也这么分析一下。不过因为 (n,i) 走不到 (n,i+1) ，所以这个题不是求无阻碍到 i=n 或 j=m 的方案数，但是到终止状态前的一个点是无阻碍的，枚举前一个点算贡献即可。

而计算 C(n,0),C(n,1),⋯,C(n,c),⋯(modp) 可以套用1.3中的做法做到 O(mlogp) 的，整体复杂度 O(k2mlogp) 。

4.5 POJ 2992 Divisors

简要题意

求 C(n,m) 的约数个数。

0≤m≤n≤431 ，保证答案不超过 (263−1)

简要题解

利用1.3中的线性筛法将 C(n,m) 表示成质数的幂次即可，答案即 ∏x is prime(cx+1) ，时间复杂度 O(n) 。

4.6 ZOJ 2116 Christopher

简要题意

给定正整数 n 和质数 p ，求 m=0,1,⋯,n 中满足 C(n,m) 是 p 的倍数的 m 个数。

1≤n≤10100,1≤p≤107

简要题解

将 n 的 p 进制按2.1的方式表示出来，可以发现 C(n,m)≢0(modp) ，则每一位都有 0≤mi≤ni ，所以不满足条件的数字个数为 ∏ki=0(ni+1) ，答案为 (n+1)−∏ki=0(ni+1) 。

4.7 HDU 5446 Unknown Treasure

简要题意

给出 n,m,k,p1,p2,⋯,pk ，保证 p1,p2,⋯,pk 是两两不同的质数，令 P=∏ki=1pi ，求 C(n,m)modP 。

1≤m≤n≤1018,1≤k≤10,1≤p1,p2,⋯,pk≤105,1≤P≤1018

简要题解

恰好符合2.2中的限定，只需要对于每个质数利用卢卡斯定理求组合数，最后再利用中国剩余定理合并即可。

4.8 CDOJ 1177 pow(a, C(n,m) ) % k

简要题意

给出 a,n,m,k ，求 aC(n,m)modk 。

1≤a≤1018,0≤n,m≤1018,1≤k≤105

简要题解

考虑 a0,a1,⋯,ak 在模 k 意义下必然存在两个相同的数字，因此可以计算出幂次的循环节 L (L≤k) （实际上 L∣φ(k) ），而 L 可以枚举得出，所以问题可以转化为 C(n,m)modL ，由于 L 不一定是square-free number，所以需要使用3.1中阶乘表示的方法。注意可能存在 a0,a1,⋯,at(modk) 不属于循环节，在  m≤n−m  时可以根据  m  的大小来确定  C(n,m)<L  的  n  范围，如果  n  真的足够小，那么可以计算实际值，否则直接计算模意义下的组合数即可。

4.9 BZOJ 4535 [Ontak2013]Kapita加强版

简要题意

给出 n,m,e ，求 C(n+m,n) 去掉所有末尾的 0 后对 10e 取模的值。

1≤n,m≤1015,1≤e≤18

简要题解

设 C(n+m,m) 的末尾有 t 个 0 ，而且 C(n+m,m)=2k1⋅v1=5k2⋅v2 ，其中 v1 与 2 互质， v2 与 5 互质。显然 t=min(k1,k2) ，因此 C(n+m,m)10t=2k1−t⋅v15t=5k2−t⋅v22t 。只要能算出 k1,v1mod2e,k2,v2mod5e ，就可以利用中国剩余定理得到答案了， 为了算出这四个数，只需要套用3.2的做法。

4.10 51nod 1387 移数字

简要题意

对于长度为 n (n≥3) 的排列定义一种操作是将第 3 个元素或最后一个元素移到排列的开头。问有多少种长度为 n 的排列能使用这种操作将其变成 1,2,⋯,n ，答案对质数 p 取模。

3≤n≤109,p=a⋅2b+1≤2⋅109 (a,b∈Z,a>0,b≥16)

简要题解

考虑操作对排列逆序对数奇偶性的变化，第一种操作的变化是偶数，第二种操作的变化与 n 的奇偶性相反，因而当 n 为偶数时答案为 n! ，当 n 为奇数时答案为 n!2 。而算阶乘的部分只需要套用3.3的做法，注意到 p=a⋅2b+1 (b≥16) ，所以FFT的复根可以利用原根的幂次来替代，从而消除误差，因此也不需要使用改良版的FFT，这里的原根是指满足 gp−1≡1(modp) 且对于任意 d>1,d|p−1 有 gp−1d≢1(modp) 的元素 g ，具体的替代方法不再赘述。

---

组合数求模的技巧暂时告一段落