看见神犇(yxs)的博客发现自己忘了(lucas),补一下
卢卡斯定理
结论:
当(p)为质数时,(C_{n}^{m}=C_{n\, mod\, p}^{m\, mod\, p}*C_{n/p}^{m/p})
证明:
引理:
((1+x)^p≡1+x^p (mod\, p))
证明:
根据费马小定理:((1+x)^p≡(1+x)≡1+x^p(mod\, p))
主体证明:
二项式定理将两侧展开
当(pi+j=m)时,(x)指数为(m),此时(i=lfloorfrac{m}{p} floor,j=m\, mod\, p)
将此时(i,j)代入右侧得到
递归求解复杂度为(O(plogp)),当(n,mge p)时只能使用(lucas)
#include<bits/stdc++.h>
using namespace std;
#define int long long
inline int read()
{
int x=0,f=1;
char ch;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
int T,n,m,p;
inline int fast(int x,int k)
{
int ret=1;
while(k)
{
if(k&1) ret=ret*x%p;
x=x*x%p;
k>>=1;
}
return ret;
}
inline int C(int n,int m)
{
if(n<m) return 0;
int sum1=1;
for(int i=n-m+1;i<=n;++i) (sum1*=i)%=p;
int sum2=1;
for(int i=1;i<=m;++i) (sum2*=i)%=p;
return sum1*fast(sum2,p-2)%p;
}
inline int lucas(int n,int m)
{
if(!m) return 1;
return lucas(n/p,m/p)*C(n%p,m%p)%p;
}
signed main()
{
T=read();
while(T--)
{
n=read(),m=read(),p=read();
printf("%lld
",lucas(n+m,m));
}
return 0;
}
拓展卢卡斯
当(p)不是质数……
数学中常见处理模数非质数套路:模数分解为质数然后(crt)
对(P)素数唯一分解(P=prodlimits_{i=1}^{k}p_i^{a_i})
求
先假设当前求的是第(k)项,即(C_{n}^{m}\, (mod\, p_k^{a_k})),后面模数用(p^k)代替
等价于(frac{n!}{(n-m)!m!}mod\, p^k)
不能直接求((n-m)!)和(n!)的逆元,因为它们不一定和(p^k)互质,可能不存在逆元
我们把左边含有(p)的因子提出来
则(frac{n!}{p^x})可以用(exgcd)求逆元
先来求(frac{n!}{p^x})中的(x)
我们知道(n)里面有(lfloorfrac{n}{p} floor)个(p)的倍数,那么我们每次统计出(lfloorfrac{n}{p} floor)个贡献递归求解
设(n!)中有(g(n))个(p)的因子,那么
边界为(g(x)=0\,(x<p))
然后考虑求(frac{n!}{p^x})的值(f(n))
由于是在模(p^k)的意义下,所以
我们把含(p)的项提出来
倒数第二个(prod)可以求一个(p^k)以内的乘积然后快速幂,最后一个(prod)直接暴力求
至于(p^{lfloorfrac{n}{p} floor})我们应该扔掉,但是((lfloorfrac{n}{p} floor)!)内可能还有答案,递归求解
复杂度为(O(plogp))
#include<bits/stdc++.h>
using namespace std;
#define int long long
inline int read()
{
int x=0,f=1;
char ch;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
int n,m,p;
inline void exgcd(int &x,int &y,int a,int b)
{
if(!b)
{
x=1,y=0;
return;
}
exgcd(y,x,b,a%b);
y-=a/b*x;
}
inline int fast(int x,int k,int p)
{
int ret=1;
while(k)
{
if(k&1) ret=ret*x%p;
x=x*x%p;
k>>=1;
}
return ret;
}
inline int inv(int a,int b)
{
int x,y;
exgcd(x,y,a,b);
return (x%b+b)%b;
}
inline int fac(int x,int a,int b)
{
if(!x) return 1;
int ret=1;
for(int i=1;i<=b;++i)
if(i%a) (ret*=i)%=b;
ret=fast(ret,x/b,b);
for(int i=1;i<=x%b;++i)
if(i%a) (ret*=i)%=b;
return ret*fac(x/a,a,b)%b;
}
inline int C(int n,int m,int a,int b)
{
int N=fac(n,a,b),M=fac(m,a,b),Z=fac(n-m,a,b);
int sum=0;
for(int i=n;i;i/=a) sum+=i/a;
for(int i=m;i;i/=a) sum-=i/a;
for(int i=n-m;i;i/=a) sum-=i/a;
return N*fast(a,sum,b)%b*inv(M,b)%b*inv(Z,b)%b;
}
inline int crt(int a,int m,int p)
{
return inv(m/p,p)*(m/p)*a;
}
inline int exlucas(int n,int m,int p)
{
int t=p,ret=0;
for(int i=2;i*i<=p;++i)
{
int k=1;
while(t%i==0) t/=i,k*=i;
if(k^1) (ret+=crt(C(n,m,i,k),p,k))%=p;
}
if(t>1) (ret+=crt(C(n,m,t,t),p,t))%=p;
return ret;
}
signed main()
{
n=read(),m=read(),p=read();
printf("%lld
",exlucas(n,m,p));
return 0;
}