概述
卢卡斯定理主要用于求解组合数取模问题。
公式
原命题见： https://zh.wikipedia.org/wiki/%E5%8D%A2%E5%8D%A1%E6%96%AF%E5%AE%9A%E7%90%86
原命题等价于：
$$\binom{m}{n}=\binom{\lfloor \dfrac{m}{p}\rfloor}{\lfloor \dfrac{n}{p}\rfloor}\times \binom{m\mod p}{n\mod p} \pmod p$$
其中 $p$ 是质数。
$\binom{m\mod p}{n\mod p}$ 中上下两项都 $< p$ ，当 $p\le 10^{7}$ 时，可以用 $O(n)$ 处理阶乘进行计算。剩下的可以不断递归进行计算。
inline int C(int n,int m)
{
if (m>n) return 0;
return fac[n]*(inv(fac[m])*inv(fac[n-m])%ha)%ha;
}
inline int lucas(int n,int m)
{
if (!m) return 1;
return lucas(n/ha,m/ha)*C(n%ha,m%ha)%ha;
}拓展卢卡斯定理
当 $p$ 很大或不是质数时，显然不能用上述方法算。
当 $p$ 不是质数，可以对其进行分解
$$p=\prod_{i=1}^n p_i^{k_i}$$
然后计算每一个 $\mod p_i^{k_i}$ 后的值，最后用中国剩余定理合并答案即可。
当 $p$ 很大时，需要对阶乘的处理进行优化。
假设当前的模数为 $p^k$ ，需要求 $n!$ 。
可以先将 $p$ 的倍数提出来，并记录个数为 $s$ 。
int s=0; //pi的个数
for (int i=n;i;i/=p) s+=i/p;$p$ 的倍数 $\div p$ 后剩下的数的乘积即为 $\lfloor \dfrac{n}{p} \rfloor !$
然后可以对剩下的数按照 $\div p^k$ 的值进行分块，长度为 $p^k$ ，计算一个块内乘积 $\mod p$ 的值就可以得到所有块的值。
还剩下 $n\mod p^k$ 个数不在块内，这些数可以单独处理。
inline int fac(int x,int pi,int pk)
{
if (!x) return 1;
int ans=1;
for (int i=2;i<=pk;i++) //处理整块
{
if (i%pi==0) continue;
ans=ans*i%pk;
}
ans=qpow(ans,x/pk,pk); //所有块的乘积
int len=x%pk;
for (int i=2;i<=len;i++) //单独处理不在块内的元素
{
if (i%pi==0) continue;
ans=ans*i%pk;
}
return ans*fac(x/pi,pi,pk)%pk; //乘上/p剩下的数
}组合数的计算公式为
$$\binom{m}{n}=\dfrac{n!}{m!\times (n-m)!}$$
显然涉及到除法，那么之前统计的 $s$ 就有用了。对于 $n!$ 就加，$m!$ 和 $(n-m)!$ 就减去 $p$ 的个数就行了。
inline int C(int n,int m,int pi,int pk)
{
int s=0; //pi的个数
for (int i=n;i;i/=pi) s+=i/pi;
for (int i=m;i;i/=pi) s-=i/pi;
for (int i=n-m;i;i/=pi) s-=i/pi;
return fac(n,pi,pk)*inv(fac(m,pi,pk),pk)%pk*inv(fac(n-m,pi,pk),pk)%pk*qpow(pi,s,pk)%pk;
}主体部分就是对 $p$ 进行分解，然后对每个因数记录 $p_i$ 为 pi，$p_i^{k_i}$ 为 pk，再对组合数进行计算，最后中国剩余定理合并答案。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
inline ll read()
{
char ch=getchar();
ll f=1,x=0;
while (ch<'0' || ch>'9')
{
if (ch=='-') f=-1;
ch=getchar();
}
while (ch>='0' && ch<='9')
{
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
ll p,zs[1005],cnt;
bool ss[1005];
inline void get_prime()
{
memset(ss,1,sizeof(ss));
ss[0]=ss[1]=0;
ll sqn=sqrt(p);
for (ll i=2;i<=sqn;i++)
{
if (ss[i]) zs[++cnt]=i;
for (ll j=1;j<=cnt && zs[j]*i<=sqn;j++)
{
ss[zs[j]*i]=0;
if (i%zs[j]==0) break;
}
}
}
ll exgcd(ll l,ll r,ll &x,ll &y)
{
if (r==0)
{
x=1; y=0;
return l;
}
ll ans=exgcd(r,l%r,y,x);
y-=l/r*x;
return ans;
}
inline ll qpow(ll x,ll y,ll ha)
{
ll ans=1;
while (y)
{
if (y&1) ans=ans*x%ha;
y>>=1;
x=x*x%ha;
}
return ans%ha;
}
inline ll inv(ll now,ll ha)
{
ll x=0,y=0;
exgcd(now,ha,x,y);
return (x%ha+ha)%ha;
}
inline ll fac(ll x,ll pi,ll pk)
{
if (!x) return 1;
ll ans=1;
for (ll i=2;i<=pk;i++) //整块
{
if (i%pi==0) continue;
ans=ans*i%pk;
}
ans=qpow(ans,x/pk,pk);
ll len=x%pk;
for (ll i=2;i<=len;i++)
{
if (i%pi==0) continue;
ans=ans*i%pk;
}
return ans*fac(x/pi,pi,pk)%pk;
}
inline ll C(ll n,ll m,ll pi,ll pk)
{
ll s=0; //pi的个数
for (ll i=n;i;i/=pi) s+=i/pi;
for (ll i=m;i;i/=pi) s-=i/pi;
for (ll i=n-m;i;i/=pi) s-=i/pi;
return fac(n,pi,pk)*inv(fac(m,pi,pk),pk)%pk*inv(fac(n-m,pi,pk),pk)%pk*qpow(pi,s,pk)%pk;
}
inline ll exlucas(ll n,ll m)
{
ll ans=0,p2=p;
for (ll i=1;i<=cnt;i++)
{
if (p2%zs[i]) continue;
ll pk=1;
while (p2%zs[i]==0) p2/=zs[i],pk*=zs[i];
ll ai=C(n,m,zs[i],pk),mi=p/pk;
ans=(ans+ai*mi%p*inv(mi,pk)%p)%p;
}
if (p2>1)
{
ll ai=C(n,m,p2,p2),mi=p/p2;
ans=(ans+ai*mi%p*inv(mi,p2)%p)%p;
}
return ans;
}
signed main()
{
ll n,m;
n=read(); m=read(); p=read();
get_prime();
return !printf("%lld",exlucas(n,m));
}