求大组合数

Posted 2020-06-29 wls001

众所周知 ,所以我们只需要求出m!和n!(m-n)!。

但是由于n和m都偏大，直接乘一定会爆掉（喜闻乐见(●‘?‘●)）。所以我们首先想到的是边乘边模，但是尝试了一组数据后我们神奇地发现答案不对。

但是我们知道另一种定理，n!%p=n!边乘边模的逆元。

所以这时候我们可以利用逆元来解决这个问题。

而逆元可以运用扩展欧几里得的方法求得。

下面我们给出欧几里得算法的证明:

首先我们设k为gcd(a,b)，则a=km，b=kn。

则a%b=a-c*b=km-c*kn=(m-cn)k

gcd(b,a%b)=gcd(kn,(m-cn)k)

由于k为a，b的最大公约数，所以n与m-cn互质，所以gcd(b,a%b)=gcd(a,b)=k。

知道了欧几里得算法，我们可以对其进行延伸得到了扩展欧几里得算法：

对于任意两个互质的a,b，总有gcd(a,b)=ax+by,扩展欧几里得算法可以用来求解

x，y。求法：

    根据欧几里得算法可知gcd(a,b)=gcd(b,a%b)。

    则bx‘+(a%b)y‘=gcd(a,b)

    将a%b=a-(a/b)*b带入得

    ay‘+b(x‘-(a/b)*y‘)=gcd(a,b)

    对于a，b而言，他们对应的x,y则分别为y‘,(x‘-(a/b)*y‘)。

而当b=0时，a*1+b*0=gcd(a,b)。

知道了这些后，此题就变得非常简单了，但是仍有需要注意的地方。

一：要防止逆元为负数的情况发生。

二：要注意因子中含有p的倍数，这样变成边模后结果就会变成0.

对于第一种情况，我们很好解决。对于第二种情况，我们可以先记录因子可以整除p几次后记录，将能整除的除去，剩下的就乘进去取模。

如果分母中的因子可以整除p的次数和与分子中的相等的话，就照原来的方法输出，否则输出0。

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

代码如下

 1 #include<iostream>
 2 #include<cstdio>
 3 #define LL long long
 4 using namespace std;
 5 LL extgcd(LL a,LL b,LL& x,LL& y)
 6 {
 7     if(b!=0)
 8     {
 9         LL d=extgcd(b,a%b,y,x);
10         y-=(a/b)*x;
11         return d;
12     }
13     else
14     {
15        x=1,y=0;
16        return a;
17     }
18 }
19 int main()
20 {
21     LL m,n,p;
22     cin>>m>>n>>p;
23     LL a=1,b=1,c=1;
24     LL ka=0,kb=0,kc=0;
25     for(int i=1;i<=m;i++)
26     {
27         int temp=i;
28          while(temp%p==0)
29           {
30             temp/=p;
31             ka++;
32           }
33           a=(a*temp)%p;
34     }
35      for(int i=1;i<=n;i++)
36      {
37         int temp=i;
38          while(temp%p==0)
39           {
40             temp/=p;
41             kb++;
42           }
43           b=(b*temp)%p;
44      }
45      for(int i=1;i<=m-n;i++)
46      {
47         int temp=i;
48          while(temp%p==0)
49           {
50             temp/=p;
51               kc++;
52           }
53           c=(c*temp)%p;
54      }
55     LL x1,x2,y;
56     extgcd(b,p,x1,y);
57     x1=(x1%p+p)%p; 
58     extgcd(c,p,x2,y);
59     x2=(x2%p+p)%p;
60     a=((a*x1)%p)*x2%p;
61     if(ka==kb+kc) cout<<a;
62     else cout<<0;
63 }

View Code