搬运自远古的洛咕博客，故文风与现在有很大不同
Lucas 卢卡斯定理
若 (p) 是质数，则对于任意整数 (1 leq m leq n) 有：
[C_n^m equiv C_{n {
m {mod}} p}^{m {
m {mod}} p} imes C_{n/p}^{m/p} ({
m {mod}} p)
]
证明比较复杂，没怎么听懂。
证明过程中有个结论可能会有用：
若 (p) 是质数，将 (a,b) 转化为 (p) 进制：
[egin{aligned}
a =a_k imes p^k +a_{k-1} imes p^{k-1}+...+a_1 imes p+ a_0\
b =b_k imes p^k +b_{k-1} imes p^{k-1}+...+b_1 imes p+ b_0\
end{aligned}
]
则：
[C_a^b equiv C_{a_k}^{b_k} imes C_{a_{k-1}}^{b_{k-1}} imes ... imes C_{a_0}^{b_0} ({
m {mod}} p)
]
模版：
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define lxl long long
using namespace std;
inline lxl pow(lxl a,lxl b,lxl p)
{
lxl ans=1%p;
while(b>0)
{
if(b&1) ans=(ans*a)%p;
a=(a*a)%p;
b>>=1;
}
return ans%p;
}
inline lxl mul(lxl a,lxl b,lxl p)
{
lxl ans=0;
while(b>0)
{
if(b&1) ans=(ans+a)%p;
a=(a+a)%p;
b>>=1;
}
return ans%p;
}
inline lxl inv(lxl x,lxl p)
{
return pow(x,p-2,p);
}
inline lxl C(lxl n,lxl m,lxl p)
{
if(n<m) return 0;
lxl up=1,down=1;
for(lxl i=n-m+1;i<=n;i++) up=mul(up,i,p);
for(lxl i=1;i<=m;i++) down=mul(down,i,p);
return mul(up,inv(down,p),p);
}
inline lxl Lucas(lxl n,lxl m,lxl p)
{
return m ? C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p :1;
}
int main()
{
//freopen("P3807.in","r",stdin);
int t;scanf("%d",&t);
lxl n,m,p;
while(t--)
{
scanf("%lld%lld%lld",&n,&m,&p);
printf("%lld
",Lucas(n+m,n,p));
}
return 0;
}
exLucas 扩展卢卡斯定理
说是扩展卢卡斯定理，但是好像和Lucas关系不大。
Lucas定理中，要求 (p) 必须为质数，那么 (p) 是任意数时怎么求呢？这里用到扩展Lucas。
设 (ans=C_n^m { m {mod}} p)。
将 (p) 分解质因数：
[p=p_1^{k_1} imes p_2^{k_2} imes ... imes p_x^{k_x}
]
则有：
[egin{cases}
ans equiv c_1 ({
m {mod}} p_1^{k_1} ) \
ans equiv c_2 ({
m {mod}} p_2^{k_2} ) \
... \
ans equiv c_x ({
m {mod}} p_x^{k_x} )
end{cases}
]
其中 (c_i=C_n^m { m {mod}} p_i^{k_i})。
也就是说，求解 (c_i) 后，再用CRT合并同余方程组即可求出 (ans)。
如何求解 (c_i) ：
[c_i=frac{n!}{m!(n-m)!} {
m {mod}} p_i^{k_i}
]
注意到 (m!,(n-m)!) 与 (p_i^{k_i}) 不一定互素，不能直接求解逆元。
考虑将 (n!,m!,(n-m)!) 中 (p_i) 因子除去，使其与 (p_i^{k_i}) 互素：
[frac{frac{n!}{p_i^x}}{frac{m!}{p_i^y} frac{(n-m)!}{p_i^z}} imes p_i^{x-y-z} {
m {mod}} p_i^{k_i}
]
其中 (x,y,z) 分别为 (n!,m!,(n-m)!) 中 (p_i) 因子的个数。
此时即可用逆元求解。
如何求解 (frac{n!}{p_i^x}) 。
对 (n!) 进行变形：
[egin{aligned}
n!&=1 imes 2 imes 3 imes ... imes n\
&=(p_i imes 2p_i imes ...)prod_{i=1,i
ot equiv 0 {
m {mod}} p_i} ^{n} i\
&=p_i^{lfloor frac {n}{p_i}
floor} imes {lfloor frac {n}{p_i}
floor}! imes prod_{i=1,i
ot equiv 0 {
m {mod}} p_i} ^{n} i
end{aligned}
]
可以发现其中 (prod_{i=1,i ot equiv 0 { m {mod}} p_i} ^{n} i) 是有循环节的，可以暴力求，例如：
[egin{aligned}
22! &equiv (1⋅2⋅4⋅5⋅7⋅8) imes (10⋅11⋅13⋅14⋅16⋅17) imes (19⋅20⋅22) imes (3⋅6⋅9⋅12⋅15⋅18⋅21)pmod {3^2}\
&equiv(1cdot 2cdot 4cdot 5cdot 7cdot 8)^2 imes (19cdot 20cdot 22) imes 3^7 imes (1cdot 2cdot 3cdot 4cdot 5cdot 6cdot 7)pmod {3^2}\
&equiv 3^7 imes 7! imes (1cdot 2cdot 4cdot 5cdot 7cdot 8)^2 imes (19cdot 20cdot 22)pmod{3^2}\
end{aligned}
]
所以：
[=p_i^{lfloor frac {n}{p_i}
floor} imes {lfloor frac {n}{p_i}
floor}! imes (prod_{i=1,i
ot equiv 0 {
m {mod}} p_i} ^{p_i^{k_i}} i)^{lfloor frac {n}{p_i^{p_k}}
floor} imes prod_{i=p_i^{k_i} imes lfloor frac {n}{p_i^{k_i}}
floor,i
ot equiv 0 {
m {mod}} p_i} ^{n} i
]
发现其中 ({lfloor frac {n}{p_i} floor}!) 和 (n!) 的求解方法是一样的，递归求解即可。
部分参考：Fading大佬的博客
模版：
#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define lxl long long
using namespace std;
inline lxl read()
{
lxl x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return x*f;
}
inline lxl fmi(lxl a,lxl b,lxl p)
{
lxl ans=1;
while(b>0)
{
if(b&1) ans=(ans*a)%p;
a=(a*a)%p;
b>>=1;
}
return ans;
}
inline lxl exgcd(lxl a,lxl b,lxl &x,lxl &y)
{
if(!b) {x=1,y=0;return a;}
lxl k=exgcd(b,a%b,x,y);
lxl z=x;x=y,y=z-a/b*y;
return k;
}
inline lxl inv(lxl a,lxl b)
{
lxl x,y;
lxl g=exgcd(a,b,x,y);
return g==1?(x+b)%b:-1;
}
inline lxl mul(lxl n,lxl pi,lxl pk)
{
if(!n) return 1;
lxl ans=1;
if(n/pk)
{
for(lxl i=1;i<=pk;i++)
if(i%pi) ans=ans*i%pk;
ans=fmi(ans,n/pk,pk);
}
for(lxl i=2;i<=n%pk;i++)
if(i%pi) ans=ans*i%pk;
return ans*mul(n/pi,pi,pk)%pk;
}
inline lxl C(lxl n,lxl m,lxl p,lxl pi,lxl pk)
{
if(n<m) return 0;
lxl a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk),k=0,ans;
for(lxl i=n;i;i/=pi) k+=i/pi;
for(lxl i=m;i;i/=pi) k-=i/pi;
for(lxl i=n-m;i;i/=pi) k-=i/pi;
ans=a*inv(b,pk)%pk*inv(c,pk)%pk*fmi(pi,k,pk)%pk;
ans=ans*(p/pk)%p*inv(p/pk,pk)%p;
return ans;
}
inline lxl exLucas(lxl n,lxl m,lxl p)
{
lxl ans=0,x=p,t=sqrt(p);
for(lxl i=2;i<=t;i++)
if(x%i==0)
{
lxl pk=1;
while(x%i==0) x/=i,pk=pk*i%p;
ans=(ans+C(n,m,p,i,pk))%p;
}
if(x>1) ans=(ans+C(n,m,p,x,x))%p;
return ans;
}
lxl n,m,p;
int main()
{
//freopen("P4720.in","r",stdin);
n=read(),m=read(),p=read();
printf("%lld
",exLucas(n,m,p));
return 0;
}