这个东西呢，需要中国剩余定理定理（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 那么第二步的a解决之后，第二步也很好办了
1 lt C(const lt n,const lt m,const lt pi,const lt pk)////C(n,m)%pk pk=pi^ki2 {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 else44 {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 完整代码
1 #include2 #include 3 #include 4 #include 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