LUCAS
(nm)≡(⌊np⌋⌊mp⌋)⋅(n%pm%p) mod  p\dbinom n m\equiv \dbinom{\lfloor\dfrac{n}{p}\rfloor}{\lfloor\dfrac{m}{p}\rfloor}·\dbinom{n\% p}{m\%p} \ \mod p(mn)≡(⌊pm⌋⌊pn⌋)⋅(m%pn%p) modp
EXLUCAS
设模数为P(P=p1k1p2k2...pnkn),piP(P=p_1^{k_1}p_2^{k_2}...p_n^{k^n}),p_iP(P=p1k1p2k2...pnkn),pi为互不相同的质数。
首先分别求出xi=(nm)mod  pikix_i=\dbinom n m\mod p_i^{k_i}xi=(mn)modpiki。
问题便转换成了求解:x≡ximod  pikix\equiv x_i\mod p_i^{k_i}x≡ximodpiki。套用CRT解决。
那么问题在于如何快速求出所有的(nm)mod  piki\dbinom n m\mod p_i^{k_i}(mn)modpiki。
当ki=1k_i=1ki=1时,套用LUCAS。
当ki>1k_i>1ki>1时,(nm)=n!m!(n−m)!\dbinom n m=\dfrac{n!}{m!(n-m)!}(mn)=m!(n−m)!n!无法对质数的次幂快速取模,考虑分别对n!,m!,(n−m)!n!,m!,(n-m)!n!,m!,(n−m)!进行拆分取模,最终答案即为n!×inv(m!,piki)×inv((n−m)!,piki)n!\times inv(m!,p_i^{k_i})\times inv((n-m)!,p_i^{k_i})n!×inv(m!,piki)×inv((n−m)!,piki)(inv(x,p)inv(x,p)inv(x,p)表示xxx在模ppp意义下的逆元,exgcd求解)。
拆分取模方式如下:
假设当前pi=3,ki=2,P=32p_i=3,k_i=2,P=3^2pi=3,ki=2,P=32,求解19!mod  P19!\mod P19!modP。
19!=(1×2×4×5×7×...×14×16×17×19)×(36)×(1×2×3×4×5×6)19!=(1\times2\times4\times5\times7\times...\times14\times16\times17\times 19)\times(3^{6})\times(1\times2\times3\times4\times5\times6)19!=(1×2×4×5×7×...×14×16×17×19)×(36)×(1×2×3×4×5×6)
对于第一部分1×2×4×5×7×8≡10×11×13×14×16×17mod  P1\times2\times 4\times 5\times 7\times 8 \equiv10\times11\times13\times14\times16\times17\mod P1×2×4×5×7×8≡10×11×13×14×16×17modP,等价于(∏i=1,i∤piP−1i)⌊nP⌋mod  P(\prod\limits_{i=1,i\nmid p_i}^{P-1}i)^{\lfloor{\frac{n}{P}}\rfloor}\mod P(i=1,i∤pi∏P−1i)⌊Pn⌋modP,剩余的数必然少于PPP个,而将右边的⌊npi⌋!\lfloor{\frac{n}{p_i}}\rfloor!⌊pin⌋!递归下去求解即可。
代码
#include<bits/stdc++.h>
#define debug(x) cerr<<#x<<":"<<x<<endl
using namespace std;
typedef long long ll;
const int N=5010;
ll n,m,p,q[N],a[N];
int tot;
namespace exlucas{
inline int fp(ll x,ll y,ll p)
{
ll re=1;
for(;y;y>>=1,x=x*x%p)
if(y&1) re=re*x%p;
return re;
}
ll exgcd(ll n,ll m,ll &x,ll &y)
{
if(!m)
{x=1;y=0;return n;}
ll re=exgcd(m,n%m,y,x);y-=(n/m)*x;
return re;
}
inline ll inv(ll n,ll p)
{
ll x,y;
exgcd(n,p,x,y);
return (x%p+p)%p;
}
inline ll C(ll n,ll m,ll p)
{
if(m>n) return 0;
ll re=1,x,y,z,i;
for(i=1;i<=m;++i){
x=n-i+1;y=(int)inv(i,p);
re=re*x%p*y%p;
}
return re;
}
ll lucas(ll n,ll m,ll p)
{
if(m==0 || n==m) return 1;
return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
ll cal(ll n,ll a,ll b,ll p)
{
if(!n) return 1;
ll i,j,re=1;
for(i=1;i<p;++i) if(i%a) re=re*i%p;
re=fp(re,n/p,p);
for(i=n/p*p+1;i<=n;++i) if(i%a) re=re*i%p;
return re*cal(n/a,a,b,p)%p;
}
inline ll multilucas(ll n,ll m,ll a,ll b,ll p)
{
ll i,s=0;
for(i=n;i;i/=a) s+=i/a;
for(i=m;i;i/=a) s-=i/a;
for(i=n-m;i;i/=a) s-=i/a;
return fp(a,s,p)*cal(n,a,b,p)%p*inv(cal(m,a,b,p),p)%p*inv(cal(n-m,a,b,p),p)%p;
}
inline ll cal(ll n,ll m,ll p)
{
ll i,j,b,c,d,ans,mod,x,y;
for(i=2;i*i<=p;++i)
if(p%i==0){
q[++tot]=1;
for(j=0;p%i==0;p/=i) q[tot]*=i,j++;
if(j==1) a[tot]=lucas(n,m,i);
else a[tot]=multilucas(n,m,i,j,q[tot]);
}
if(p>1){q[++tot]=p;a[tot]=lucas(n,m,p);}
ans=a[1];mod=q[1];
for(i=2;i<=tot;++i){
b=exgcd(mod,q[i],x,y);
c=a[i]-ans;
if(c%b) return -1LL;
d=q[i]/b;x=((c/b)%d*x%d+d)%d;
ans=ans+mod*x;
mod=mod/b*q[i];
}
return ans;
}
}
int main(){
scanf("%lld%lld%lld",&n,&m,&p);
printf("%lld\n",exlucas::cal(n,m,p));
return 0;
}