卢卡斯定理与拓展卢卡斯

看见神犇(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))

主体证明：

[(1+x)^n≡(1+x)^{lfloorfrac{n}{p} floor p}*(1+x)^{n\, mod\, p}(mod\, p) ]

[(1+x)^n≡(1+x^p)^{lfloorfrac{n}{p} floor}*(1+x)^{n\, mod\, p}(mod\, p) ]

二项式定理将两侧展开

[sumlimits_{i=0}^{n}C_{n}^{i}x^i≡sumlimits_{i=0}^{lfloorfrac{n}{p} floor}C_{lfloorfrac{n}{p} floor}^{i}x^{pi}*sumlimits_{i=0}^{n\, mod\, p}C_{n\, mod\, p}^{j}x^j(mod\, p) ]

当(pi+j=m)时，(x)指数为(m)，此时(i=lfloorfrac{m}{p} floor,j=m\, mod\, p)

将此时(i,j)代入右侧得到

[C_{n}^{m}=C_{lfloorfrac{n}{p} floor}^{lfloorfrac{m}{p} floor}*C_{n\, mod\, p}^{m\, mod\, p}(mod\, p) ]

递归求解复杂度为(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})

求

[egin{cases}C_{n}^{m}\, (mod p_1^{a_1})\C_{n}^{m}\, (mod p_2^{a_2})\……\C_{n}^{m}\, (mod p_n^{a_n})end{cases} ]

先假设当前求的是第(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{frac{n!}{p^x}}{frac{m!}{p^y}frac{(n-m)!}{p^z}}p^{x-y-z}mod\, p^k ]

则(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(n)=lfloorfrac{n}{p} floor+g(lfloorfrac{n}{p} floor) ]

边界为(g(x)=0\,(x<p))

然后考虑求(frac{n!}{p^x})的值(f(n))

由于是在模(p^k)的意义下，所以

[n!=(prodlimits_{i=1}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1}^{n\, mod\, p^k}i ]

我们把含(p)的项提出来

[n!=prodlimits_{i=1}^{lfloorfrac{n}{p} floor}(pi)(prodlimits_{i=1,i\,mod\,p e 0}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1,i\,mod\,p e 0}^{n\, mod\, p^k}i ]

[n!=p^{lfloorfrac{n}{p} floor}(lfloorfrac{n}{p} floor)!(prodlimits_{i=1,i\,mod\,p e 0}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1,i\,mod\,p e 0}^{n\, mod\, p^k}i ]

倒数第二个(prod)可以求一个(p^k)以内的乘积然后快速幂，最后一个(prod)直接暴力求

至于(p^{lfloorfrac{n}{p} floor})我们应该扔掉，但是((lfloorfrac{n}{p} floor)!)内可能还有答案，递归求解

[f(n)=f(lfloorfrac{n}{p} floor)(prodlimits_{i=1,i\,mod\,p e 0}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1,i\,mod\,p e 0}^{n\, mod\, p^k}i ]

复杂度为(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;
}