组合数计算（从1000到1e9的组合数各类求法）

最新推荐文章于 2022-05-11 22:04:37 发布

wizardAEI

最新推荐文章于 2022-05-11 22:04:37 发布

阅读量2.2k

点赞数 6

分类专栏：积累算法篇文章标签：算法

本文链接：https://blog.csdn.net/qq_37360324/article/details/108868134

版权

积累算法篇专栏收录该内容

8 篇文章 0 订阅

订阅专栏

组合数涉及的问题很多，求解过程中涉及的定理判定也有很多，期间断断续续参考了很多大佬们的博客，最终进行了总结。在此我使用用问题引入算法的方法来一步一步了解组合数。

简单的组合数问题（运用快速幂及递归实现杨辉三角）

组合数求解

像关于组合数求和问题如∑ C_nⁱ（0<=i<=n）可利用快速幂，此时∑ C_nⁱ = 2^n

此为二项式定理推导出的结果练习题可见洛谷P3414
若求较小的组合数如C(n,m)（m<=n<=1000 保守估计1000，实际可能可以更大些，但不会超数量级）可直接利用杨辉三角的法则来计算组合数C(n,m)=C(n-1,m-1)+C(n-1，m)。以下为杨辉三角的简单介绍和理解：

杨辉三角简单介绍

性质

1.每个数等于它上方两数之和–>最重要的一个性质

2.每行数字左右对称，由1开始逐渐变大–>由第5个性质可知

3.第n行的数字有n项

4.第n行的m个数可表示为C(n-1,m-1),即为从n-1个不同元素中取m-1个元素的组合数

5.第n行的第m个数和第n-m+1个数相等,为组合数性质之一,即C(n,m)==C(n,n-m+1)

6.每个数字等于上一行的左右两个数字之和。可用此性质写出整个杨辉三角。即第n+1行的第i个数等于第n行的第i-1个数和第i个数之和，这也是组合数的性质之一。即 C(n+1,i)=C(n,i)+C(n,i-1)。

7.(a+b)n的展开式中的各项系数依次对应杨辉三角的第(n+1)行中的每一项–>二项式定理

杨辉三角还有其他性质（例如斐波那契数列与其的关系）在此不做说明

由第6个性质可知上面第二类问题C(n,m)=C(n-1,m-1)+C(n-1，m)的来由

但是，利用杨辉三角递归，只可以处理n,m较小的组合数。当n，m很大时，我们如何求解呢?

较复杂的组合数问题(n,m较大时)

逆元在组合数中的应用

对欧拉-费马小定理的理解

在了解逆元之前，需要先了解欧拉-费马小定理，并有此推导逆元（如果已经了解了欧拉-费马小定理，可直接跳到标题：通过逆元计算组合数）

欧拉定理：若正整数 a , n 互质（公因数只有1），则 a^φ(n)≡1(mod n) 其中 φ(n) 是欧拉函数，此函数代表（1~n) 与 n 互质的数的个数。

简单证明:

令 x₁ x₂ x₃…x_φ(n) 为1-n中与n互质的数（互质即gcd(m,n)=0)

我们先来研究 ax₁ ax₂ ax₃… ax_φ(n) 这一队数列

**其必有任意两个元素mod n 不相同 **

反证：若存在ax1 ≡ ax2（mod n），则 n |（ ax1 - ax2 ）即n可以整除( ax1 - ax2 )

而（ax1-ax2）/n ==（x1-x2）/n * a 或 a/n * (x1-x2) ；

第一种除法由于（x1-x2）< n ∴一定除不开

第二种除法由于 a和 n互质 ∴也一定除不开。综上不可能存在ax1 ≡ ax2（mod n），得证。

且∵ a与n互质所以 ax_i 与n也互质所以根据欧几里得算法（辗转相除法）有gcd（axi，n）1 所以：gcd（ax_i，n） gcd（n，ax_i%n）== 1

又∵ax₁%n, ax₂%n, ax₃%n… ax_φ(n)%n 均小于 n，且ax₁%n, ax₂%n, ax₃%n… ax_φ(n)%n又都与n互质

∴ ax₁%n, ax₂%n, ax₃%n… ax_φ(n)%n必为x₁ x₂ x₃…x_φ(n) 以一定顺序排列而成的。

那么必有：ax₁%n, ax₂%n, ax₃%n… ax_φ(n)%n = x₁ x₂ x₃…x_φ(n)

自然就有 ax₁%n, ax₂%n, ax₃%n… ax_φ(n)%n ≡ x₁ x₂ x₃…x_φ(n) （mod n）

即 ax₁ ax₂ ax₃ … ax_φ(n) ≡ x₁ x₂ x₃…x_φ(n) （mod n），把x₁ x₂ x₃…x_φ(n) 移到左边（因为肯定不会出现小数，所以可以直接移）

∴有 **(a^φ(n)-1)x₁ x₂ x₃…x_φ(n) ≡ 0 (mod n) **

又∵x₁ x₂ x₃…x_φ(n) 均与 n互质

∴自然有(a^φ(n)-1) ≡ 0 (mod n)

∴ 证得 a^φ(n)≡1(mod n)

接下来，我们利用欧拉定理，简单证明费马小定理

费马小定理：对于质数p，任意整数，均满足：a^p≡a（mod p）

简单证明：

由于 a^p≡a（mod p）∴有a^p-1≡1 (mod p) 这和 a^φ(n)≡1(mod n)类似，所以只需证明p-1φ(n)即可，且因为p是一个质数，易得与其互质的数为1，2，3…p-1,共p-1个所以有p-1φ(n)，得证

(注意：若有a是p的倍数，则一定有a^p≡a≡0（mod p），不必证明)

（拓展）欧拉定理的推论：若正整数a，n互质，那么对于任意正整数b，有a^b≡a^{b mod φ（n）}（mod n）

简单证明：若a^b ≡ a^{b mod φ（n）}（mod n)则有a^{b-b mod φ（n）}≡ 1（mod n)，因为φ（n）| b-b mod φ（n）所以令 p=（b-b mod φ（n））/φ（n） ∴即证(a^p)^φ（n）≡1 (mod n) 易得a^p和 n是互质的所以得证

以上引用，总结自：https://www.cnblogs.com/zylAK/p/9569668.html （这篇博客里还有一道例题，想具体了解欧拉定理推论可以去看一下）

通过逆元计算组合数

首先了解逆元的定义：

逆元：对于a和p（a和p互素），若a*b%p≡1，则称b为a的逆元。(mod p的情况下)

首先，我们弄清，为什么求组合数需要用到逆元呢？

此时我们就需要使用上面👆说过的费马小定理。

根据欧拉-费马小定理可知 a^φ§ ≡ 1 (mod p) 又因为当p为质数时，满足欧拉定理，且φ§=p-1，所以有a^p-2 × a≡1(mod p),所以有a^p-2 是a的逆元，此时，我们只需要把a^p-2 求出来，即可求出逆元。然后须用到快速幂算法（快速幂相关内容可在我的博客主页搜索：快速幂算法）

finally，让我们来总结下如何求：
$C_n^m\mod p (p为质数)$

求出n!%p，m!%p，(n-m)!%p

(注1：阶乘过程中，由同模运算的乘法原则，可在每一步阶乘中都mod上p)

(注2：若求同n下的多个m，即求多组组合数，可以定义F[i],为i！mod p的值，那么就可以把多组的m！%p和(n-m)！%p 都存在F数组中)
求出m!^p-2%p和(n-m)!^p-2%p （即对m!%p，(n-m)!%p 快速幂）
求出(n!%p×m!^p-2%p*(n-m)!^p-2%p)%p 即为最终答案。

以上引用，总结自：https://www.cnblogs.com/liziran/p/6804803.html

例题：牛客挑战赛37 B

简单概述：输入一组合数C_n^m,若与所给数相同，输出Yes!否则输出No!

思路：使用逆元求解，但这里面没有给出模数，我们可以采用随机质数的方法，打一个随机质数表，通过多次取余与所给数取余比较，来判断是否相同。注意到所给数p可能很大，在保存时，使用字符串保存，再分部取余。

（这是一个AC代码的质数表，可以适当采用： { 88897 ,58787 ,125453 , 7853 , 78569 ,145267 ,25867 ,105899 ,8803,4591 }）

👇下给出通过逆元求组合数的做法

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int mo[10]={88897,58787},mt,nt;
ll mp=0;
char p[400010];
void mop(int mod)
{
    mp=0;
    int len=strlen(p);
    for(int i=0;i<len-1;i++)
    mp=(mp+(p[i]-'0')%mod)*10%mod;mp=(mp+p[len-1]-'0')%mod;
}
ll fastpow(ll baser,ll power,int mod)
{
    ll ans=1;
    while(power)
    {
        if(power&1)
        ans=ans*baser%mod;
        baser=baser*baser%mod,power>>=1;
    }
    return ans%mod;
}
ll ie(ll n,ll m,int mod)//inverse element 的缩写
{
    ll nj=1,mj=1,nmj=1;
    for(int i=2;i<=m;i++)
        mj=mj*i%mod;
    for(int i=2;i<=n-m;i++)
        nmj=nmj*i%mod;
        nj=mj;
    for(int i=m+1;i<=n;i++)
        nj=nj*i%mod;
    return (nj*fastpow(mj,mod-2,mod)%mod*fastpow(nmj,mod-2,mod)%mod)%mod;
}
int main()
{

    cin>>nt>>mt>>p;
    for(int i=0;i<2;i++)
    {
        mop(mo[i]);//求mp
        if(ie(nt,mt,mo[i])!=mp){cout<<"No!";return 0;}
    }
    cout<<"Yes!";
}

那么，我们会得到以下结果：

分析为什么会超时呢？

我们查看题目描述，发现1≤n≤1e9；

在进行阶乘运算时，由于其复杂度为O(n)，那么可知，当n接近于1e9时，复杂度会过大

那么剩余的60％，我们如何去求组合数？那么我们就需要了解Lucas定理👇

Lucas定理在组合数中的应用

(注意以下的除法均为向下取整)

*Lucas定理的定义：C(n,m)%p=C(n/p,m/p)C(n%p,m%p)%p

那么我们只要不断对C(n/p,m/p)调用Lucas定理就可以通过递归来实现对组合数计算的简化,其复杂度为O(logⁿ)

卢卡斯定理证明

(~~网上找了很多证明都不太懂，去了b站听，豁然开朗😂~~)

首先我们令：n=a×p+b m=c×p+d，也就是a=n/p b=n%p c=m/p d=m%p

那么Lucas定理可以写成：
$C_n^m\equiv C_a^b*C_c^d (\mod p)$
我们先把这种写法放一边，来看这两个式子：

(1+x)^p ≡ 1+x （mod p）以及 1+x^p ≡ 1+x （mod p）

这两个式子我们都可以用欧拉—费马小定理来证明，所以有(1+x)^p ≡ 1+x^p

那么这个同模有啥用呢我们可以利用它来得到下面这个式子👇(n=a×p+b m=c×p+d)

(1+x)ⁿ = (1+x)^a×p+b

= (1+x)^p×a × (1+x)^b

≡ (1+x^p)^a × (1+x)^b (mod p)

接下来，我们利用二项式定理，可知道在式子左边(1+x)ⁿ 中的x^m 项为：
$C_n^mx^m$
而式子右端 (1+x^p)^a × (1+x)^b 的x^m 项为
$C_a^c C_b^d x^{p*c}x^d=C_a^c C_b^d x^{p*c+d}=C_a^c C_b^d x^{m}$
又∵式子左端和右端同余

∴有
$C_n^m≡C_a^c C_b^d (\mod p)$
得证。

🐱‍👤进阶：Lucas定理的真正面目（看懂上面的递归就可以解决问题了）

我们由上面可以知道：C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p

然后C(n/p,m/p)继续套用Lucas定理，那么，我们真正在一直乘的其实是后面的C(n%p,m%p)，且这个这里面的两个余数，是经过某些次连除p后得到的"n"和"m"再取余

那么我们将n，m写成p进制数，则有

n=(n₁n₂n₃…n_i)_p m=(m₁m₂m₃…m_j)_p

我们令两个p进制数右对齐，m若位数和n不同，用0在左边补上，保证都为K位数。

那么这p进制下的n,m相对应的每一位数，其实就是一个个n%p,m%p！

那么我们就可以得知Lucas定理的另一种写法(原始写法)：
$C_n^m≡\prod_{i=1}^{K} {C_{n_i}^{m_i}} \mod p$

具体使用方法

Lucas(n,m)=Lucas(n/p,m/p)*ie(n%p,m%p) 其中ie为逆元函数

边界问题：lucas和ie函数 m==0时返回1，ie函数 n<m时返回 0；

所以，我们继续完善以上的例题

例题AC代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int mo[10]={88897,58787,125453,7853,78569,145267,25867,105899,8803,4591},mt,nt;
ll mp=0;
char p[400010];
void mop(int mod)
{
    mp=0;
    int len=strlen(p);
    for(int i=0;i<len-1;i++)
    mp=(mp+(p[i]-'0')%mod)*10%mod;mp=(mp+p[len-1]-'0')%mod;
}
ll fastpow(ll baser,ll power,int mod)
{
    ll ans=1;
    while(power)
    {
        if(power&1)
        ans=ans*baser%mod;
        baser=baser*baser%mod,power>>=1;
    }
    return ans%mod;
}
ll ie(ll n,ll m,int mod)//inverse element 的缩写
{
    if(n<m) return 0;if(m==0) return 1;
    ll nj=1,mj=1,nmj=1;
    for(int i=2;i<=m;i++)
        mj=mj*i%mod;
    for(int i=2;i<=n-m;i++)
        nmj=nmj*i%mod;
        nj=mj;
    for(int i=m+1;i<=n;i++)
        nj=nj*i%mod;
    return (nj*fastpow(mj,mod-2,mod)%mod*fastpow(nmj,mod-2,mod)%mod)%mod;
}
ll lucas(ll n,ll m,int mod)
{
    if(m==0) return 1;
    return (lucas(n/mod,m/mod,mod)*ie(n%mod,m%mod,mod))%mod;
}
int main()
{

    cin>>nt>>mt>>p;
    for(int i=0;i<10;i++)
    {
        mop(mo[i]);//求mp
        if(lucas(nt,mt,mo[i])!=mp){cout<<"No!";return 0;}
    }
    cout<<"Yes!";
}

以上，我们都是在讨论模数为质数时的情况，那么当模数不为质数时，我们就要用到扩展Lucas定理👇

扩展Lucas定理(当模数不为质数时)：

中国剩余定理：

为了了解并证明扩展lucas定理，我们先了解什么是中国剩余定理：

1585043571712

其证明可以在知乎中得到答案：https://zhuanlan.zhihu.com/p/44591114 (写得很好)

拓展Lucas定理具体使用：

当C_n^m mod p 其中p不为质数时，那么对p进行质因子分解，有
$p=\prod_i p_i^{k_i}$
其中k_i表示含有的第i个质因子的个数，比如12=2²×3¹。那么我们易得p_i^k_i 之间是互素的

接下来，我们只要分别求出
$a_i≡C_n^m\mod p_i^{k_i}$
对应上方👆中国剩余定理我们可知，x就是C_n^m ，且：
$C_n^m≡ \sum_i (a_i\times \frac p {p_i^{k_i}}\times[ (\frac p {p_i^{k_i}} )^- ]_ {p_i^{k_i}}) \mod p$
其中
$(\frac p {p_i^{k_i}} )^- ]_ {p_i^{k_i}})$
表示分式的逆元 mod p_i^k_i (~~看不懂就去看看上面的知乎链接，比对来理解~~)

那么我们只剩下最终一个问题：如何求a_i ?

以下图引用自:https://www.bilibili.com/video/BV1xb411x74F?from=search&seid=7989126052916540468

由于
$a_i≡C_n^m\mod p_i^{k_i}$
则问题可以转换为以下图片的求解：

其中k1 k2 k3 分别为 n! m! (n-m)! 的阶乘中所含的因子p的个数

→得知一般规律后带入，（(m!/p^k₂)%p^k)^{φ(p^k)-1} ,（((n-m)!/p^k₃)%p^k)^{φ(p^k)-1} 就是逆元

我们以n=19,p=3,k=2为例：

n!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19

=(1∗2∗4∗5∗7∗8）∗（10∗11∗13∗14∗16∗17）∗19∗3⁶∗(1∗2∗3∗4∗5∗6)

将其分为3部分：第一部分是p的幂的部分，也就是3⁶即p^⌊n/p⌋；第二部分是一个新的阶乘，也就是6!即⌊n/p⌋!，并且可以继续递归直到其中不含有质因子p;第三部分是除了前两部分之外剩下的数

我们注意到p^a全部由第一部分和第二部分提供，且a为n!中所有的质因子p的个数；

所以(n!/p^a)%p^k 就是第二部分递归所得的最终的值乘上第三部分 mod p^k（程序中顺便记录出a的值好最后再乘回来）

那么第三部分有什么特点呢？

注意到(1∗2∗4∗5∗7∗8）≡（10∗11∗13∗14∗16∗17）（mod p^k）所以直接求第一组再快速幂就可以了。幂值为组数，为⌊n/p^k⌋。

且19<=p^k 可以暴力直接乘上

第三部分怎么来的呢？

首先第一组为(1∗2∗4∗5∗7∗8)是数i(1<=i<p^k)且与p互质的数，且和其他组数加起来总组数为⌊n/p^k⌋

这个19 就是最终剩下的没有形成完整一组的数其个数易得为从1到n%p^k中与p互质的个数（记作s）；且其相乘的模数就是第一组的前 s个数相乘模p^k的值

即

那么
$C_n^m \mod p$
向回推导

（(m!/p^k₂)%p^k)^{φ(p^k)-1} ,（((n-m)!/p^k₃)%p^k)^{φ(p^k)-1} 此逆元就可以求出来了,然后根据“通过逆元计算组合数”，就可以算出最终答案。

p不为质数的组合数求解，到此其中的所有问题都已解决,可以求出最终答案🐱‍👤。

由于菜鸡本菜还没遇到相关类型的题，无法测试扩展Lucas算法的代码正确性，所以没有贴出来，有兴趣者可以在其他地方找到相关代码。

终言：洋洋洒洒五千字，说不尽组合事。欲穷千里目，还需更上一层楼。

安利一个markdown编辑器：Typora，真的好用。~~其实也是偶然发现的( •̀ ω •́ )y~~

我的个人博客网站aei.today ~~但是GitHub托管好像出了问题，可能会找不到ip地址~~

wizardAEI

关注

6
点赞
踩
11

收藏

觉得还不错? 一键收藏
打赏
1
评论
组合数计算（从1000到1e9的组合数各类求法）

组合数涉及的问题很多，求解过程中涉及的定理判定也有很多，期间断断续续参考了很多大佬们的博客，最终进行了总结。在此我使用用问题引入算法的方法来一步一步了解组合数。简单的组合数问题（运用快速幂及递归实现杨辉三角）组合数求解像关于组合数求和问题如∑ Cni（0<=i<=n）可利用快速幂，此时∑ Cni = 2^n此为二项式定理推导出的结果练习题可见洛谷P3414若求较小的组合数如C(n,m)（m<=n<=1000 保守估计1000，实际可能可以更大些，但不会超数量级
复制链接

扫一扫