題目傳送門

求組合數的時候，如果模數p是質數，可以用盧卡斯定理解決。

但是盧卡斯定理僅僅適用於p是質數的情況。

當p不是質數的時候，我們就需要用擴充套件盧卡斯求解。

實際上，擴充套件盧卡斯=快速冪+快速乘+exgcd求逆元+質因數分解+crt合併答案+求階乘，跟盧卡斯定理沒什麼關係......

如果把模數p分解成p1^k1*p2^k2*...*px^kx的形式，那麼我們可以求出c(n,m)分別模每個pi^ki的結果，再用中國剩餘定理合併即可。

每個pi^ki一定是互質的，所以用樸素crt就行。

根據組合數的定義，c(n,m)=(n!) / (m!*(n-m)!) ，所以我們只要能想辦法求出階乘，就能再利用exgcd求出逆元，進而求出組合數。

接下來唯一的問題就是怎麼快速求出 x! 取模 pi^ki 的結果。

考慮如下的經典樣例（據說來自popoqqq）：(19!)%(3^2)

19!=1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19

先把其中的3的倍數提出來，因為求組合數的時候分子分母能約掉。

19!=(1*2*4*5*7*8)*(10*11*13*14*16*17)*(19)*(3*6*9*12*15*18)=(1*2*4*5*7*8)*(1*2*4*5*7*8)*(3*3*3*3*3*3)*(1*2*3*4*5*6)=(1*2*4*5*7*8)^2*19*(3^6)*(1*2*3*4*5*6)。

後面的6!部分可以遞迴求解，遞迴終點為0!=1。

3^6最後計算組合數的時候再處理。

那幾個(1*2*4*5*7*8)顯然是迴圈的，迴圈節長度小於pi^ki，可以暴力計算。

顯然一共有(x/(pi^ki))個迴圈節，套個快速冪即可。

剩下的部分，即19，長度等於x%(pi^ki)，也小於pi^ki，也可以暴力計算。

至此我們求出了階乘。

求組合數的時候，考慮pi的倍數的影響。

分子分母分別計數相加減。

最後用crt合併即可。

 1 #include<cstdio>
 2 typedef long long ll;
 3 
 4 ll n,m,p;
 5
 
 
 6 ll ksm(ll b,ll tp,ll mod)
 7 {
 8     ll ret=1;
 9     while(tp)
10     {
11         if(tp&1)ret=ret*b%mod;
12         b=b*b%mod;
13         tp>>=1;
14     }
15     return ret;
16 }
17 
18 ll mul(ll a,ll b,ll mod)
19 {
20     ll ret=0;
21     while(b)
22     {
23         if(b&1)ret=(ret+a)%mod;
24         a=(a+a)%mod;
25         b>>=1;
26     }
27     return ret;
28 }
29 
30 ll exgcd(ll a,ll b,ll &x,ll &y)
31 {
32     if(!b)
33     {
34         x=1;y=0;
35         return a;
36     }
37     ll t=exgcd(b,a%b,y,x);
38     y-=a/b*x;
39 }
40 
41 ll inv(ll x,ll mod)
42 {
43     ll a,b;
44     exgcd(x,mod,a,b);
45     return (a%mod+mod)%mod;
46 }
47 
48 ll fac(ll x,ll pi,ll pk)
49 {
50     if(!x)return 1;
51     ll ans=1;
52     for(ll i=2;i<=pk;i++)
53         if(i%pi)ans=ans*i%pk;
54     ans=ksm(ans,x/pk,pk);
55     for(ll i=2;i<=x%pk;i++)
56         if(i%pi)ans=ans*i%pk;
57     return ans*fac(x/pi,pi,pk)%pk;
58 }
59 
60 ll c(ll cn,ll cm,ll pi,ll pk)
61 {
62     if(cm>cn)return 0;
63     ll up=fac(cn,pi,pk),d1=fac(cm,pi,pk),d2=fac(cn-cm,pi,pk);
64     ll cnt=0;
65     for(ll i=cn;i;i/=pi)cnt+=i/pi;
66     for(ll i=cm;i;i/=pi)cnt-=i/pi;
67     for(ll i=cn-cm;i;i/=pi)cnt-=i/pi;
68     return up*inv(d1,pk)%pk*inv(d2,pk)%pk*ksm(pi,cnt,pk)%pk;
69 }
70 
71 ll crt(ll a,ll pk)
72 {
73     return a*inv(p/pk,pk)%p*(p/pk)%p;
74 }
75 
76 int main()
77 {
78     scanf("%lld%lld%lld",&n,&m,&p);
79     ll tp=p,ans=0;
80     for(ll i=2;i*i<=p;i++)
81     {
82         if(tp%i)continue;
83         ll pk=1;
84         while(!(tp%i))tp/=i,pk*=i;
85         ans=(ans+crt(c(n,m,i,pk),pk))%p;
86     }
87     if(tp>1)ans=(ans+crt(c(n,m,tp,tp),tp))%p;
88     printf("%lld",(ans%p+p)%p);
89     return 0;
90 }