233333,本来是想自己写的,可写到一半发现有一个不会的地方,又去现学的,蓝后,百度到一篇非常好的博客,自认为简单易懂,博主提供的模板也好记,不如直接转载qwq
首先给出Lucas(卢卡斯)定理:
有非负整数A、B,和素数p,A、B写成p进制为:A=a[n]a[n-1]...a[0],B=b[n]b[n-1]...b[0]。
则组合数C(A,B)与C(a[n],b[n])×C(a[n-1],b[n-1])×...×C(a[0],b[0]) mod p同余。
即:Lucas(n,m,p)=C(n%p,m%p)×Lucas(n/p,m/p,p) ,特别的,Lucas(x,0,p)=1。
其实说白了,Lucas定理就是求组合数C(n,m)mod p(p是素数)的值,即(( n! /(n-m)!)/m! )mod p,而我们
又知道(a/b)mod p=a*b^(p-2)mod p(这里用到了一点逆元和费马小定理的知识),这样我们就可以在计算阶乘的过程中对p取模,不会造成溢出。(需要注意的是Lucas定理处理的p的范围大致为10^5数量级)
看下面一道题:HDU3037
题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=3037
题目大意:m个种子要放在n颗树上,问有多少种方法,结果对素数p取模。
分析:m个种子可以分成两部分:放在树上的和不放在树上的,我们可以假想出第n+1颗树,认为那些没放在树上的种子放在这颗假想树上,这样,问题就变为了m个种子放在n颗树上有多少种方案数了。等价于方程X1+X2+...+Xn+1=m有多少组非负整数解,即(X1+1)+(X2+1)+...+(Xn+1+1)=m+n+1有多少组正整数解。挡板原理:(m+n+1)个1,分成n+1部分的方案数==>C(n+m,n)。到这儿就很明显了,果断Lucas。
代码:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<iostream> 2 #include<cstdio> 3 #define ll long long 4 using namespace std; 5 const int maxn=150010; 6 int t; 7 ll n,m,p; 8 ll jc[maxn]; 9 void work(ll mod) 10 { 11 jc[0]=1; 12 for(int i=1;i<=p;++i) 13 jc[i]=jc[i-1]*i%mod; 14 } 15 ll qpow(ll a,ll b) 16 { 17 ll ans=1; 18 while(b) 19 { 20 if(b&1)ans=(ans*a%p)%p; 21 a=(a*a)%p; 22 b>>=1; 23 } 24 return ans; 25 } 26 ll C(ll n,ll m) 27 { 28 if(m>n)return 0; 29 return jc[n]*qpow(jc[m]*jc[n-m]%p,p-2)%p;//定理,a%p的逆元等于a^(p-2)%p 30 } 31 ll lucas(ll n,ll m) 32 { 33 if(m==0)return 1;//c(n,0)=1; 34 else return (C(n%p,m%p)*lucas(n/p,m/p)%p)%p; 35 } 36 int main() 37 { 38 scanf("%d",&t); 39 while(t--) 40 { 41 scanf("%lld%lld%lld",&n,&m,&p); 42 work(p); 43 printf("%lld ",lucas(n+m,m)); 44 } 45 return 0; 46 }
对于p比较大的情况,不能对阶乘预先处理,需要单独处理。
题目链接:http://acm.fzu.edu.cn/problem.php?pid=2020
题目大意:求C(n,m)mod p。n,p是10^9数量级的。
代码:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <iostream> 2 #include <string.h> 3 #include <stdio.h> 4 5 using namespace std; 6 typedef long long LL; 7 8 LL n,m,p; 9 10 LL quick_mod(LL a, LL b) 11 { 12 LL ans = 1; 13 a %= p; 14 while(b) 15 { 16 if(b & 1) 17 { 18 ans = ans * a % p; 19 b--; 20 } 21 b >>= 1; 22 a = a * a % p; 23 } 24 return ans; 25 } 26 27 LL C(LL n, LL m) 28 { 29 if(m > n) return 0; 30 LL ans = 1; 31 for(int i=1; i<=m; i++) 32 { 33 LL a = (n + i - m) % p; 34 LL b = i % p; 35 ans = ans * (a * quick_mod(b, p-2) % p) % p; 36 } 37 return ans; 38 } 39 40 LL Lucas(LL n, LL m) 41 { 42 if(m == 0) return 1; 43 return C(n % p, m % p) * Lucas(n / p, m / p) % p; 44 } 45 46 int main() 47 { 48 int T; 49 scanf("%d", &T); 50 while(T--) 51 { 52 scanf("%I64d%I64d%I64d", &n, &m, &p); 53 printf("%I64d ", Lucas(n,m)); 54 } 55 return 0; 56 }