Lucas定理是用来求 c(n,m) mod p,p为素数的值。
- 表达式:
- C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p
- 当我们遇到求一个N,M很大的组合数的时候,递推法就显得很耗时了,对于1e9那么大的数据求N!%P,无论是空间还是时间都不会允许。
- 于是引申出lucas定理,利用这个表达式可以将数量级降低好几个,从而减小时间和空间的开销。
- 一般来说可以用lucas定理解决的问题都是N,M很大,但质数P相对来说在1e5左右,不会太大,我们利用迭代渐渐缩小N,M的值,
- 将C(N%P,M%P)累乘在答案上求解,下面来思考怎么编写代码。
这个定理目前我不会证明,只是知道内容,惭愧。
当N,M不为0时且组合数合法我们可以继续迭代,当算出来的N%P<M%P时表示结果为0此时直接返回0即可。
对于一个较小的C(n,m)=n!/(m!*(n-m)!) ,我们就可以根据这个式子得到 C(n,m)%P=f(n)%P*(f(m)*f(n-m))-1%P 其中f表示阶乘
由于P是一个质数,由费马小定理可以得到 (f(m)*f(n-m))-1=mod_pow(f(m)*f(n-m),P-2,P),我们可以先打表出f[P]之内的阶乘对P取余然后直接调用即可。
下面是主要代码:
1 //lucas定理 2 LL quick_mod(LL a, LL b, LL c) 3 { 4 LL ans = 1; 5 while(b) 6 { 7 if(b & 1) 8 ans = (ans*a)%c; 9 b>>=1; 10 a = (a*a)%c; 11 } 12 return ans; 13 } 14 LL fac[MAXN_P]; 15 void Get_Fac(LL p)///m! 16 { 17 fac[0] = 1; 18 for(int i=1; i<=p; i++) 19 fac[i] = (fac[i-1]*i) % p; 20 } 21 LL Lucas(LL n, LL m, LL p) //C(n,m)%p 22 { 23 LL ans = 1; 24 while(n && m) 25 { 26 LL a = n % p; 27 LL b = m % p; 28 if(a < b) 29 return 0; 30 ans = ( (ans*fac[a]%p) * (quick_mod(fac[b]*fac[a-b]%p,p-2,p)) ) % p; 31 n /= p; 32 m /= p; 33 }
下面以51nod 1120为例 , 这是题目链接
是一个经典的求卡特兰数的例题,但是N达到了1e9的规模,如果使用O(N)打表,f[n]=(f[n-1]*(4*n-2))/(n+1)得话时间不允许。
还有另一种递推式子是 f[n]=C(N*2,N)/(N+1) ,所以答案就是二倍的第(N-1)个卡特兰数对10007取余,可以用lucas求解,注意有分母的存在要求一下逆元,
由于N很大,也不能用打表法求逆元,可以用拓展欧几里得或者快速幂(由于这里P是质数且lucas就要用到快速幂,所以快速幂很方便不必要重写exgcd)
参考代码:
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define LL long long 4 #define MAXN 100005 5 void gcd(LL a,LL b,LL &d,LL &x,LL &y) 6 { 7 if(!b) {d=a;x=1;y=0;} 8 else {gcd(b,a%b,d,y,x);y-=x*(a/b);} 9 } 10 LL Inv(LL a,LL n) 11 { 12 LL d,x,y; 13 gcd(a,n,d,x,y); 14 return d==1?(x+n)%n:-1; 15 } 16 17 LL quick_mod(LL a, LL b, LL c) 18 { 19 LL ans = 1; 20 while(b) 21 { 22 if(b & 1) 23 ans = (ans*a)%c; 24 b>>=1; 25 a = (a*a)%c; 26 } 27 return ans; 28 } 29 LL fac[MAXN]; 30 void Get_Fac(LL p)///m! 31 { 32 fac[0] = 1; 33 for(int i=1; i<=p; i++) 34 fac[i] = (fac[i-1]*i) % p; 35 } 36 LL Lucas(LL n, LL m, LL p) 37 { 38 LL ans = 1; 39 while(n && m) 40 { 41 LL a = n % p; 42 LL b = m % p; 43 if(a < b) 44 return 0; 45 ans = ( (ans*fac[a]%p) * (quick_mod(fac[b]*fac[a-b]%p,p-2,p)) ) % p; 46 n /= p; 47 m /= p; 48 } 49 return ans; 50 } 51 int main() 52 { 53 LL n, m, p=10007; 54 Get_Fac(p); 55 cin>>n;n--; 56 cout<<Lucas(n*2,n,p)*2*quick_mod(n+1,p-2,p)/*Inv(n+1,p)*/%p<<endl; 57 return 0; 58 }