扩展lucas定理
前置知识:
- lucas定理
- 中国剩余定理
介绍
当正整数n,mn,mn,m很大,且质数ppp较小的时候,要求CnmC_n^mCnm对ppp取模后的值,可以用lucas定理。
但如果ppp不是质数,那该怎么办呢?如果mmm较小,则可以用扩展lucas定理。
第一步:中国剩余定理
设p=p1r1p2r2⋯pkrkp=p_1^{r_1}p_2^{r_2}\cdots p_k^{r_k}p=p1r1p2r2⋯pkrk,其中pip_ipi为质数。我们可以先求出Cnm%p1r1,Cnm%p2r2,…,Cnm%pkrkC_n^m\%p_1^{r_1},C_n^m\%p_2^{r_2},\dots,C_n^m\%p_k^{r_k}Cnm%p1r1,Cnm%p2r2,…,Cnm%pkrk的值a1,a2,…,aka_1,a_2,\dots,a_ka1,a2,…,ak。
我们把CnmC_n^mCnm看作未知数xxx,可以得到以下方程组:
{x≡a1(modp1r1)x≡a2(modp2r2)x≡a3(modp3r3)......x≡an(modpkrk)\left\{ \begin{matrix} x\equiv a_1\pmod{p_1^{r_1}}\\ x\equiv a_2\pmod{p_2^{r_2}}\\ x\equiv a_3\pmod{p_3^{r_3}}\\ ......\\ x\equiv a_n\pmod{p_k^{r_k}} \end{matrix} \right. ⎩⎨⎧x≡a1(modp1r1)x≡a2(modp2r2)x≡a3(modp3r3)......x≡an(modpkrk)
利用中国剩余定理,我们可以求出xxx,它是以ppp为周期出现的无穷多个解。而在[0,p)[0,p)[0,p)这个周期的解,就是Cnm%pC_n^m\%pCnm%p后的值。
那么a1,a2…,aka_1,a_2\dots,a_ka1,a2…,ak怎么求呢?
第二步:组合数模质数的幂
由第一步可得
a=Cnmmodpra=C_n^m\bmod p^ra=Cnmmodpr
因为Cnm=n!m!(n−m)!C_n^m=\dfrac{n!}{m!(n-m)!}Cnm=m!(n−m)!n!,我们若要求m!m!m!和(n−m)!(n-m)!(n−m)!关于prp^rpr的逆元,则要把其中所有的质因子ppp提出来,再乘回去即可。
Cnm=n!m!(n−m)!=n!pxm!py×(n−m)!pz×px−y−zC_n^m=\dfrac{n!}{m!(n-m)!}=\dfrac{\frac{n!}{p^x}}{\frac{m!}{p^y}\times \frac{(n-m)!}{p^z}}\times p^{x-y-z}Cnm=m!(n−m)!n!=pym!×pz(n−m)!pxn!×px−y−z
其中x,y,zx,y,zx,y,z分别是n!,m!,(n−m)!n!,m!,(n-m)!n!,m!,(n−m)!中质因子ppp的次数。此时m!py×(n−m)!pz\dfrac{m!}{p^y}\times \dfrac{(n-m)!}{p^z}pym!×pz(n−m)!与prp^rpr互质,可以直接求逆元。因为CnmC_n^mCnm为整数,所以x−y−z≥0x-y-z\geq 0x−y−z≥0,px−y−zp^{x-y-z}px−y−z可以用快速幂来求。
第三步:阶乘除去质因子后模质数幂
接下来的问题就是计算以下式子
n!ptmodpk\dfrac{n!}{p^t}\bmod p^kptn!modpk
我们呢先考虑如如何计算n!modpkn!\bmod p^kn!modpk。举个例子:n=22,p=3,k=2n=22,p=3,k=2n=22,p=3,k=2
22!=1×2×3×4×5×6×7×8×9×10×11×12×13×14×15×16×17×18×19×20×21×2222!=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10\times 11\times 12\times 13\times 14\times 15\times 16\times 17\times 18\times 19\times 20\times 21\times 2222!=1×2×3×4×5×6×7×8×9×10×11×12×13×14×15×16×17×18×19×20×21×22
把其中333的倍数提出来,得到
22!=(3×6×9×12×15×18×21)×(1×2×4×5×7×8×10×11×13×14×16×17×19×20×22)22!=(3\times 6\times 9\times 12\times 15\times 18\times 21)\times (1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22)22!=(3×6×9×12×15×18×21)×(1×2×4×5×7×8×10×11×13×14×16×17×19×20×22)
=37×(1×2×3×4×5×6×7)×(1×2×4×5×7×8×10×11×13×14×16×17×19×20×22)\qquad =3^7\times (1\times 2\times 3\times 4\times 5\times 6\times 7)\times (1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22)=37×(1×2×3×4×5×6×7)×(1×2×4×5×7×8×10×11×13×14×16×17×19×20×22)
其中373^737即为pkp^kpk,就是需要被提出的部分。
对于7!7!7!,即为⌊np⌋!\lfloor \dfrac np\rfloor!⌊pn⌋!,可以递归来求。
对于后面的部分,我们发现
1×2×4×5×7×8≡10×11×13×14×16×17(modpk)1\times 2\times 4\times 5\times 7\times 8\equiv 10\times 11\times 13\times 14\times 16\times 17\pmod{p^k}1×2×4×5×7×8≡10×11×13×14×16×17(modpk)
我们发现1×2×4×5×7×81\times 2\times 4\times 5\times 7\times 81×2×4×5×7×8在整个式子中会出现⌊npk⌋\lfloor\dfrac{n}{p^k}\rfloor⌊pkn⌋次,因此,我们可以先计算在pkp^kpk以内的部分,然后再求其⌊npk⌋\lfloor\dfrac{n}{p^k}\rfloor⌊pkn⌋次幂。不要忘了乘上最后多出的一部分。
1×2×4×5×7×8×10×11×13×14×16×17×19×20×22≡(1×2×4×5×7×8)3×19×20×22(modpk)1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22\equiv (1\times 2\times 4\times 5\times 7\times 8)^3\times 19\times 20\times 22\pmod{p^k}1×2×4×5×7×8×10×11×13×14×16×17×19×20×22≡(1×2×4×5×7×8)3×19×20×22(modpk)
也就是说,对于以下式子
=37×(1×2×3×4×5×6×7)×(1×2×4×5×7×8×10×11×13×14×16×17×19×20×22)\qquad =3^7\times (1\times 2\times 3\times 4\times 5\times 6\times 7)\times (1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22)=37×(1×2×3×4×5×6×7)×(1×2×4×5×7×8×10×11×13×14×16×17×19×20×22)
373^737是要提出的,不用计算。第二部分可以递归计算。第三部分可以O(pk)O(p^k)O(pk)得出。
总结
扩展lucas定理与lucas定理在实现上并没有太大关联,只是解决的问题比较类似。扩展lucas定理的时间复杂度大概在O(p+log2n)O(p+\log^2 n)O(p+log2n)。当然,这是最坏的时间复杂度,一般的时间复杂度远远低于此。如果ppp的质因子比较多且都比较小,则时间复杂度会降低很多。
例题
P4720 【模板】扩展卢卡斯定理
code
#include<bits/stdc++.h>
using namespace std;
int tot=0;
long long mod,x,y,ans=0,a[105],r[105];
long long mi(long long t,long long v){if(v==0) return 1;long long re=mi(t,v/2);re=re*re%mod;if(v&1) re=re*t%mod;return re;
}
void exgcd(long long c,long long d){if(d==0){x=1;y=0;return;}exgcd(d,c%d);long long t=x;x=y;y=t-c/d*y;
}
long long gt(long long v,long long p,long long q){if(!v) return 1;long long re=1;for(int i=1;i<=q;i++){if(i%p) re=re*i%q;}re=mi(re,v/q)%q;for(int i=1;i<=v%q;i++){if(i%p) re=re*i%q;}return re*gt(v/p,p,q)%q;
}//第三步
long long C(long long v1,long long v2,long long p,long long q){if(v1<v2) return 0;long long f1=gt(v1,p,q),f2=gt(v2,p,q),f3=gt(v1-v2,p,q),vt=0;for(long long i=p;i<=v1;i*=p) vt+=v1/i;for(long long i=p;i<=v2;i*=p) vt-=v2/i;for(long long i=p;i<=v1-v2;i*=p) vt-=(v1-v2)/i;return mi(p,vt)%q*f1%q*(mi(f2,q-q/p-1)%q)%q*(mi(f3,q-q/p-1)%q)%q;
}//第二步
int main()
{long long n,m,v;scanf("%lld%lld%lld",&n,&m,&mod);v=mod;for(int i=2;i*i<=v;i++){if(v%i==0){r[++tot]=1;while(v%i==0){r[tot]*=i;v/=i;}a[tot]=C(n,m,i,r[tot]);}}if(v>1){r[++tot]=v;a[tot]=C(n,m,v,v);}v=mod;for(int i=1;i<=tot;i++){exgcd(v/r[i],r[i]);x=(x%r[i]+r[i])%r[i];ans=(ans+v/r[i]*a[i]*x%v)%v;}//第一步printf("%lld",ans);return 0;
}