扩展卢卡斯
upd:2020.10.15 重写
一.什么是扩展Lucas
解决组合数对非质数取模。
二.算法流程
由于P可能为合数，先将P分解质因数。分解为(P=p_1^{k_1}*p_2^{k_2}...p_n^{k_n})对于任意两个(p_i^{k_i}),是互质的，并且有几个方程。形如
[C_n^mequiv C_n^m\%p_i^{k_i}(mod;p_i^{k_i})
]
可用中国剩余定理（CRT）将各个方程合并（不同质因子之间肯定互质，模数乘起来就会变作P）。此处挂一个算法主体与CRT
ll CRT(ll r,ll pi){return r*Inv(mod/pi,pi)%mod*(mod/pi)%mod;}
ll Exlucas(ll a,ll b){
ll tmp=mod,res=0;
for(ll i=2;i*i<=mod;++i)
if(tmp%i==0){
ll pk=1;
while(!(tmp%i))tmp/=i,pk*=i;
(res+=CRT(C(a,b,i,pk),pk))%=mod;
}
if(tmp>1)(res+=CRT(C(a,b,tmp,tmp),tmp))%=mod;
return res%mod;
}
(mod为模数)
三.阶乘的本质
[n!=p^{⌊frac{n}{p}⌋}*⌊frac{n}{p}⌋!* ∏_{p∤i}^{n}i
]
其中 p 可取任何值。（即使大于n也可以，想想为什么）
很好理解，就只是选了一个数 p ，将 (1-n) 中的 p 的倍数全部提取出来，然后再将系数提出来，分成系数相乘和 p 的幂次，剩下不是p的倍数的原样处理。
举个例子，
[egin{align}
&1*2*3*4*5*6*7*8*9\
=&(3*6*9)*(1*2*4*5*7*8)\
=&(1*3*2*3*3*3)*(1*2*4*5*7*8)\
=&(1*2*3)*3^3*(1*2*4*5*7*8)
end{align}
]
这个柿子具有可以递归的性质，这很重要。至于有什么具体的用处请见下文。
四.求解组合数模单个质数乘方的值
先把组合数的式子展开。
[frac{n!}{(n-m)!m!}\%p_i^{k_i}
]
分母目前无法用逆元处理，因为并不与 ({p_i}^{k_i}) 互质。为了达到互质的目的，将所有的 (p_i) 因子提取出来，其中 (f(n,p_i)) 为 (n!) 中因子 (p_i) 的个数，于是。
[frac{n!}{(n-m)!m!}equivfrac{frac{n!}{p_i^{f(n,p_i)}}}{frac{m!}{p_i^{f(m,p_i)}}*frac{(n-m)!}{p_i^{f(n-m,p_i)}}}*p_i^{f(n,p_i)-f(m,p_i)-f(n-m,p_i)} pmod{p_i^{k_i}}
]
已经将多余的因子影响去掉了，可以逆元。此时需要一个求解(x!~\%~p^k)的函数，并在其中不乘上 (p) 的倍数。
挂一个C的写法。（其中jc指求解(x!~\%~p^k)的函数）
ll C(ll a,ll b,ll pi,ll pk){
ll k=f(a,pi)-f(b,pi)-f(a-b,pi);
if(!KSM(pi,k,pk))return 0;
return jc(a,pi,pk)*Inv(jc(b,pi,pk),pk)%pk*Inv(jc(a-b,pi,pk),pk)%pk*KSM(pi,k,pk)%pk;
}
至于 f 既可以用递归，也可以非递归。递归形式为
[f(n,p)=f(lfloorfrac{n}{p}
floor,p)+lfloorfrac{n}{p}
floor
]
（可以参见三大点进行理解，非递归放在第五大点中）
四.写出阶乘函数（去掉 p 的影响）
[n!=p^{⌊frac{n}{p}⌋}*⌊frac{n}{p}⌋!* ∏_{p∤i}^{}n_i
]
首先柿子很明显有递归性质。含 (p) 的第一项舍去（不计p的影响），第二项递归。考虑第三项。
因为是对于 ({p_i}^{k_i}) 进行取模，可以对于要乘的一系列连续的数进行划分，长度为 ({p_i}^{k_i}) 。用柿子来解释（第一个结合律展开一下消掉含有(p^k)的柿子）。
[egin{align}
ecauseprodlimits_{i=u*{p}^{k}+1,p∤i}^{(u+1)*{p}^{k}}iequivprodlimits_{i=u*{p}^{k}+1,p∤i}^{(u+1)*{p}^{k}}i+p^k pmod{p^k}
\ hereforeprodlimits_{i=u*{p}^{k}+1,p∤i}^{(u+1)*{p}^{k}}iequivprodlimits_{i=(u+1)*{p}^{k}+1,p∤i}^{(u+2)*{p}^{k}}i pmod{p^k}
end{align}
]
划分之后只取一段然后快速幂即可。注意处理最后凑不齐一段的数。
ll jc(ll a,ll pi,ll pk){
if(!a)return 1;
ll res=1;
m_for(i,2,pk)if(i%pi)res=(res*i)%pk;
res=KSM(res,a/pk,pk);
m_for(i,2,a%pk)if(i%pi)res=(res*i)%pk;
return res*jc(a/pi,pi,pk)%pk;
}
五.Code
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define m_for(i,a,b)for(int i=a;i<=b;++i)
const int N=50010;
ll mod;
ll n,k;
ll f(ll a,ll pi){return a?f(a/pi,pi)+(a/pi):0;}
/*
ll f(ll a,ll pi){
ll res=0;
while(a)res+=a/pi,a/=pi;
return res;
}
*/
ll Exgcd(ll a,ll b,ll &x,ll &y){
if(!b)x=1,y=0;
else Exgcd(b,a%b,y,x),y-=a/b*x;
}
ll Inv(int a,int p){
ll x,y;Exgcd(a,p,x,y);
return (x%p+p)%p;
}
inline ll KSM(ll a,ll b,ll p){
ll res=1;
for(;b;b>>=1,(a*=a)%=p)if(b&1)(res*=a)%p;
return res%p;
}
ll jc(ll a,ll pi,ll pk){
if(!a)return 1;
ll res=1;
m_for(i,2,pk)if(i%pi)res=(res*i)%pk;
res=KSM(res,a/pk,pk);
m_for(i,2,a%pk)if(i%pi)res=(res*i)%pk;
return res*jc(a/pi,pi,pk)%pk;
}
ll C(ll a,ll b,ll pi,ll pk){
ll k=f(a,pi)-f(b,pi)-f(a-b,pi);
if(!KSM(pi,k,pk))return 0;
return jc(a,pi,pk)*Inv(jc(b,pi,pk),pk)%pk*Inv(jc(a-b,pi,pk),pk)%pk*KSM(pi,k,pk)%pk;
}
ll CRT(ll r,ll pi){return r*Inv(mod/pi,pi)%mod*(mod/pi)%mod;}
ll Exlucas(ll a,ll b){
ll tmp=mod,res=0;
for(ll i=2;i*i<=mod;++i)
if(tmp%i==0){
ll pk=1;
while(!(tmp%i))tmp/=i,pk*=i;
(res+=CRT(C(a,b,i,pk),pk))%=mod;
}
if(tmp>1)(res+=CRT(C(a,b,tmp,tmp),tmp))%=mod;
return res%mod;
}
int main(){
scanf("%lld%lld%lld",&n,&k,&mod);
printf("%lld",Exlucas(n,k));
return 0;
}