【EXLUCAS模板】【拓展卢卡斯详解】【组合数高级篇】LuoGu P4720
这个东西呢，需要中国剩余定理定理（crt）和拓展欧几里得（exgcd）组合数还有逆元的知识
我们要求的呢，就是这个东西\[C_n^m\% p\]
这里P不一定是质数
1.将P质因数分解
\[p = \prod\limits_{\rm{i}}^r {P_{}^{{c_i}}} \]
对下面这个进行中国剩余定理合并，便是答案
\[\begin{array}{l}
\begin{array}{*{20}{l}}
{{x_1} \equiv C_n^m(\,\bmod \,p_1^{{c_1}})}\\
{{x_2} \equiv C_n^m(\,\bmod \,p_2^{{c_2}})}\\
{{x_3} \equiv C_n^m(\,\bmod \,p_3^{{c_3}})}
\end{array}\\
\cdot \cdot \cdot \\
{x_r} \equiv C_n^m(\,\bmod \,p_r^{{c_r}})
\end{array}\]
2.求$C_n^m\% {p^k}$
这个该怎么求呢，但是很明显这个求不了逆元，因为分母的两个数和模数不互质，那么只要把所有的模数因子全部干掉就好了，这样就能求逆元了
也就是这样\[C_n^m\% {p^k} = \frac{{\frac{{n!}}{{{p^{{a_1}}}}}}}{{\frac{{m!}}{{{p^{{a_2}}}}}*\frac{{(n - m)!}}{{{p^{{a_3}}}}}}}*{p^{{a_1} - {a_2} - {a_3}}}\% {p^k}\]
但是这里的$a1a2a3$是不知道的，没事我们看下一步就知道了
3.求$n!\% {p^k}$
我们来举个栗子
\[\begin{array}{l}
22!\% {3^2} = (1*2*3*4*5*6*7*8*9*10*11*12*13*14*14*16*17*18*19*20*21*22)\% {3^2}\\
= (3*6*9*12*15*18*21)*(1*2*4*5*7*8)*(10*11*13*14*16*17)*(19*20*22)\% {3^2}\\
= {3^7}*7!*{(1*2*4*5*7*8)^2}\% {3^2}
\end{array}\]
我们推广一下\[n!\% {p^k} = {p^{n/p}}*(n/p)!*{(\prod\limits_1^{ < = {p^k}} {nu{m_i}} )^{n/{p^k}}}*(\prod\limits_1^{p\% k} {nu{m_i}} )\% {p^k}\]
看到$(n/p)!$，这个很明显可以递归来解决（注意边界是$n=1||n=2$时候返回值是1）
看到那个${p^{n/p}}$这不就是第二步里没解决的a值吗，对，就是$n/p$；
看那个${(\prod\limits_1^{ < = {p^k}} {nu{m_i}} )^{n/{p^k}}}$,这个只要写一个循环，计算出一个循环节，然后再快速幂即可；
看那个$(\prod\limits_1^{p\% k} {nu{m_i}} )$这个是冗余，写个循环就可以计算
但是$n/p$很明显是递归的，这个应该怎么解决，写个循环每次都除，再累加就好了
就是这样
1 lt fac(const lt n,const lt pi,const lt pk)//n! % pi^ki pk=pi^ki 2 { 3 if(n==1||n==0)return 1; 4 lt ans=1; 5 for(rt i=2;i<pk;i++) 6 if(i%pi) 7 ans=ans*i%pk; 8 ans=ksm(ans,n/pk,pk); 9 for(rt i=2;i<=n%pk;i++) 10 if(i%pi)ans=ans*i%pk; 11 return ksc(ans,fac(n/pi,pi,pk),pk); 12 }
那么第二步的a解决之后，第二步也很好办了
1 lt C(const lt n,const lt m,const lt pi,const lt pk)////C(n,m)%pk pk=pi^ki 2 { 3 lt up=fac(n,pi,pk),dl=fac(m,pi,pk),dr=fac(n-m,pi,pk),kk=0; 4 for(rt i=n;i;i/=pi)kk+=i/pi; 5 for(rt i=m;i;i/=pi)kk-=i/pi; 6 for(rt i=n-m;i;i/=pi)kk-=i/pi; 7 return up*inv(dl,pk)%pk*inv(dr,pk)%pk*ksm(pi,kk,pk)%pk; 8 }
然后我进行了一些基本函数丧心病狂的封装23333
1 namespace basic 2 { 3 inline lt read() 4 { 5 rt x = 0; char zf = 1; char ch; 6 while (ch != '-' && !isdigit(ch)) ch = getchar(); 7 if (ch == '-') zf = -1, ch = getchar(); 8 while (isdigit(ch)) x = x * 10 + ch - '0', ch = getchar(); return x * zf; 9 } 10 void write(lt x) 11 { 12 if (x<0)putchar('-'),x=-x; 13 if (x>9) write(x/10); 14 putchar(x%10+'0'); 15 } 16 inline void writeln(const lt x){write(x);putchar('\n');} 17 lt gcd(lt a,lt b){return b?gcd(b,a%b):a;} 18 lt ksc(lt a,lt b,lt mod) 19 { 20 lt fina=0,kk=1; 21 if(a<0)a=-a,kk=-kk;if(b<0)b=-b,kk=-kk; 22 while(b) 23 { 24 if(b%2)fina=(fina+a)%mod; 25 a=(a+a)%mod,b>>=1; 26 } 27 return fina%mod*kk; 28 } 29 lt ksm(lt a,lt k,lt mod) 30 { 31 if(!k)return 1; 32 lt fina=1; 33 while(k) 34 { 35 if(k%2)fina=ksc(fina,a,mod); 36 a=a*a%mod,k>>=1; 37 } 38 return fina%mod; 39 } 40 void ex_gcd(lt a,lt b,lt &x,lt &y) 41 { 42 if(!b)x=1,y=0; 43 else 44 { 45 ex_gcd(b,a%b,x,y);lt tt=x;x=y,y=tt-a/b*x; 46 } 47 } 48 lt exgcd(lt a,lt b,lt c)//ax=c(mod b) 49 { 50 lt gc=gcd(a,b),x=0,y=0;; 51 if(c%gc)return -1; 52 a/=gc,b/=gc,c/=gc; 53 ex_gcd(a,b,x,y); 54 return (x*c%b+b)%b; 55 } 56 inline lt inv(lt a,lt p) 57 { 58 if(!a)return 0; 59 return exgcd(a,p,1); 60 }//ax=1(mod p) 61 }
然后lucas部分
1 namespace lucas 2 { 3 using namespace basic; 4 lt fac(const lt n,const lt pi,const lt pk)//n! % pi^ki pk=pi^ki 5 { 6 if(n==1||n==0)return 1; 7 lt ans=1; 8 for(rt i=2;i<pk;i++) 9 if(i%pi) 10 ans=ans*i%pk; 11 ans=ksm(ans,n/pk,pk); 12 for(rt i=2;i<=n%pk;i++) 13 if(i%pi)ans=ans*i%pk; 14 return ksc(ans,fac(n/pi,pi,pk),pk); 15 } 16 lt C(const lt n,const lt m,const lt pi,const lt pk)////C(n,m)%pk pk=pi^ki 17 { 18 lt up=fac(n,pi,pk),dl=fac(m,pi,pk),dr=fac(n-m,pi,pk),kk=0; 19 for(rt i=n;i;i/=pi)kk+=i/pi; 20 for(rt i=m;i;i/=pi)kk-=i/pi; 21 for(rt i=n-m;i;i/=pi)kk-=i/pi; 22 return up*inv(dl,pk)%pk*inv(dr,pk)%pk*ksm(pi,kk,pk)%pk; 23 } 24 void exlucas() 25 { 26 n=read(),m=read(),P=read(); 27 lt ans=0,p=P; 28 for(rt i=2;i<=p;i++) 29 if(p%i==0) 30 { 31 lt pi=i,pk=1; 32 while(!(p%i))p/=i,pk*=i; 33 ans=(ans+C(n,m,pi,pk)*(P/pk)%P*inv(P/pk,pk)%P)%P; 34 } 35 writeln(ans%P); 36 } 37 }
完整代码
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cctype> 5 typedef long long lt; 6 #define rt register lt 7 using namespace std; 8 lt n,m,p,P; 9 namespace basic 10 { 11 inline lt read() 12 { 13 rt x = 0; char zf = 1; char ch; 14 while (ch != '-' && !isdigit(ch)) ch = getchar(); 15 if (ch == '-') zf = -1, ch = getchar(); 16 while (isdigit(ch)) x = x * 10 + ch - '0', ch = getchar(); return x * zf; 17 } 18 void write(lt x) 19 { 20 if (x<0)putchar('-'),x=-x; 21 if (x>9) write(x/10); 22 putchar(x%10+'0'); 23 } 24 inline void writeln(const lt x){write(x);putchar('\n');} 25 lt gcd(lt a,lt b){return b?gcd(b,a%b):a;} 26 lt ksc(lt a,lt b,lt mod) 27 { 28 lt fina=0,kk=1; 29 if(a<0)a=-a,kk=-kk;if(b<0)b=-b,kk=-kk; 30 while(b) 31 { 32 if(b%2)fina=(fina+a)%mod; 33 a=(a+a)%mod,b>>=1; 34 } 35 return fina%mod*kk; 36 } 37 lt ksm(lt a,lt k,lt mod) 38 { 39 if(!k)return 1; 40 lt fina=1; 41 while(k) 42 { 43 if(k%2)fina=ksc(fina,a,mod); 44 a=a*a%mod,k>>=1; 45 } 46 return fina%mod; 47 } 48 void ex_gcd(lt a,lt b,lt &x,lt &y) 49 { 50 if(!b)x=1,y=0; 51 else 52 { 53 ex_gcd(b,a%b,x,y);lt tt=x;x=y,y=tt-a/b*x; 54 } 55 } 56 lt exgcd(lt a,lt b,lt c)//ax=c(mod b) 57 { 58 lt gc=gcd(a,b),x=0,y=0;; 59 if(c%gc)return -1; 60 a/=gc,b/=gc,c/=gc; 61 ex_gcd(a,b,x,y); 62 return (x*c%b+b)%b; 63 } 64 inline lt inv(lt a,lt p) 65 { 66 if(!a)return 0; 67 return exgcd(a,p,1); 68 }//ax=1(mod p) 69 } 70 namespace lucas 71 { 72 using namespace basic; 73 lt fac(const lt n,const lt pi,const lt pk)//n! % pi^ki pk=pi^ki 74 { 75 if(n==1||n==0)return 1; 76 lt ans=1; 77 for(rt i=2;i<pk;i++) 78 if(i%pi) 79 ans=ans*i%pk; 80 ans=ksm(ans,n/pk,pk); 81 for(rt i=2;i<=n%pk;i++) 82 if(i%pi)ans=ans*i%pk; 83 return ksc(ans,fac(n/pi,pi,pk),pk); 84 } 85 lt C(const lt n,const lt m,const lt pi,const lt pk)////C(n,m)%pk pk=pi^ki 86 { 87 lt up=fac(n,pi,pk),dl=fac(m,pi,pk),dr=fac(n-m,pi,pk),kk=0; 88 for(rt i=n;i;i/=pi)kk+=i/pi; 89 for(rt i=m;i;i/=pi)kk-=i/pi; 90 for(rt i=n-m;i;i/=pi)kk-=i/pi; 91 return up*inv(dl,pk)%pk*inv(dr,pk)%pk*ksm(pi,kk,pk)%pk; 92 } 93 void exlucas() 94 { 95 n=read(),m=read(),P=read(); 96 lt ans=0,p=P; 97 for(rt i=2;i<=p;i++) 98 if(p%i==0) 99 { 100 lt pi=i,pk=1; 101 while(!(p%i))p/=i,pk*=i; 102 ans=(ans+C(n,m,pi,pk)*(P/pk)%P*inv(P/pk,pk)%P)%P; 103 } 104 writeln(ans%P); 105 } 106 } 107 int main(){lucas::exlucas();return 0;}
