组合数取模（转载）
本文转自：http://blog.csdn.net/skywalkert/article/details/52553048
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-01-01
1. 枚举法
这一部分的内容主要是利用基本定义
来进行分析的。
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)，这是因为
如果能提前O(m)预处理1,2,⋯,m的逆元，则上述问题可以做到O(m)，这是因为
注意，当m不为质数时，(mmodx)不一定与m互质，所以上述做法需要有m是质数的条件。
1.2 预处理
由于
，所以能够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)完成组合数的质因数分解。具体做法如下：
- 用筛法计算出正整数x(x≤n)的一个质因数prx。
- 维护一个数组c1,c2,⋯,cn使得C(n,m)=∏ni=1ici，初始化每个位置为0，然后使c1,c2,⋯,cn各加一，c1,c2,⋯,cm各减一，c1,c2,⋯,cn−m各减一。
- 降序枚举x(2≤x≤n)，如果x不为质数，则将xcx拆分为prxcx⋅(xprx)cx，具体操作是修改cx,cprx,cxprx。
- 此时cx不为0则必然x为质数，而且x为质数时cx非负，即为组合数的质因数分解。
当然，如果只有少量数字没有逆元，还可以考虑更改模数。例如，只有x没有逆元，但是要除以一个x，则可以计算原式在模x⋅p意义下的值，再除以x，而如果有多个数字时，也可将模数扩大它们的乘积倍。但一般情况下，这种方法的效率不高。
2. 定理法
这一部分的内容将涉及初等数论中两个常用的定理，利用定理将复杂的问题简单化，从而便于分析。
2.1 卢卡斯定理 (Lucas’s theorem)
在数论中，卢卡斯定理将组合数C(n,m)在模质数p意义下的值表示成一系列组合数C(ni,mi)的乘积，其中ni,mi分别表示n,m在p进制下的第i位数字，换句话说，令n=∑ki=0ni⋅pi,m=∑ki=0mi⋅pi，则有
其中，如果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)
中国剩余定理适用于化简模线性方程组，具体来说，如果有
，其中gcd(mi,mj)=1(i≠j)，那么利用中国剩余定理可以得到
与原方程组等价，其中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意义下进行运算。
考虑这样一个事实：如果已知[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≥e3>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 更复杂的情况
可以发现，当e>1时3.2的做法比3.3的做法更好，但也不难发现，3.3中提到的分块技巧是可以与3.2中的分治技巧进行结合的，如果对于更加一般化的做法感兴趣，可以参考Andrew Granville在1997年发表的文章Arithmetic Properties of Binomial Coefficients I: Binomial coefficients modulo prime powers。
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，则有
，对应的方案数是
。对于最后的求和部分，先考虑i=n的情况，则0≤j≤m的情况为
，那么也是枚举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=prx(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，具体的替代方法不再赘述。
