Loading

总结「卢卡斯定理」

搬运自远古的洛咕博客,故文风与现在有很大不同


Lucas 卢卡斯定理

\(p\) 是质数,则对于任意整数 \(1 \leq m \leq n\) 有:

\[C_n^m \equiv C_{n {\rm \ {mod}} \ p}^{m {\rm \ {mod}} \ p} \ \times \ C_{n/p}^{m/p} \ ({\rm {mod}} \ p) \]

证明比较复杂,没怎么听懂。

证明过程中有个结论可能会有用:

\(p\) 是质数,将 \(a,b\) 转化为 \(p\) 进制:

\[\begin{aligned} a =a_k \times p^k +a_{k-1} \times p^{k-1}+...+a_1 \times p+ a_0\\ b =b_k \times p^k +b_{k-1} \times p^{k-1}+...+b_1 \times p+ b_0\\ \end{aligned} \]

则:

\[C_a^b \equiv C_{a_k}^{b_k} \times C_{a_{k-1}}^{b_{k-1}} \times ... \times C_{a_0}^{b_0} \ ({\rm {mod}} \ p) \]


模版:

P3807 【模板】卢卡斯定理

#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\n",Lucas(n+m,n,p));
	}
	return 0;
}


exLucas 扩展卢卡斯定理

说是扩展卢卡斯定理,但是好像和Lucas关系不大。

Lucas定理中,要求 \(p\) 必须为质数,那么 \(p\) 是任意数时怎么求呢?这里用到扩展Lucas。

\(ans=C_n^m \ {\rm {mod}} \ p\)

\(p\) 分解质因数:

\[p=p_1^{k_1} \times p_2^{k_2} \times ... \times p_x^{k_x} \]

则有:

\[\begin{cases} ans \equiv c_1 ({\rm {mod}} \ p_1^{k_1} ) \\ ans \equiv c_2 ({\rm {mod}} \ p_2^{k_2} ) \\ ... \\ ans \equiv c_x ({\rm {mod}} \ p_x^{k_x} ) \end{cases} \]

其中 \(c_i=C_n^m \ {\rm {mod}} \ p_i^{k_i}\)

也就是说,求解 \(c_i\) 后,再用CRT合并同余方程组即可求出 \(ans\)

如何求解 \(c_i\)

\[c_i=\frac{n!}{m!(n-m)!} \ {\rm {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}} \times p_i^{x-y-z} \ {\rm {mod}} \ p_i^{k_i} \]

其中 \(x,y,z\) 分别为 \(n!,m!,(n-m)!\)\(p_i\) 因子的个数。

此时即可用逆元求解。


如何求解 \(\frac{n!}{p_i^x}\)

\(n!\) 进行变形:

\[\begin{aligned} n!&=1 \times 2 \times 3 \times ... \times n\\ &=(p_i \times 2p_i \times ...)\prod_{i=1,i \not \equiv 0 \ {\rm {mod}} \ p_i} ^{n} i\\ &=p_i^{\lfloor \frac {n}{p_i} \rfloor} \times {\lfloor \frac {n}{p_i} \rfloor}! \times \prod_{i=1,i \not \equiv 0 \ {\rm {mod}} \ p_i} ^{n} i \end{aligned} \]

可以发现其中 \(\prod_{i=1,i \not \equiv 0 \ {\rm {mod}} \ p_i} ^{n} i\) 是有循环节的,可以暴力求,例如:

\[\begin{aligned} 22! &\equiv (1⋅2⋅4⋅5⋅7⋅8)\times (10⋅11⋅13⋅14⋅16⋅17)\times (19⋅20⋅22)\times (3⋅6⋅9⋅12⋅15⋅18⋅21)\pmod {3^2}\\ &\equiv(1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)^2\times (19\cdot 20\cdot 22)\times 3^7\times (1\cdot 2\cdot 3\cdot 4\cdot 5\cdot 6\cdot 7)\pmod {3^2}\\ &\equiv 3^7\times 7!\times (1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)^2\times (19\cdot 20\cdot 22)\pmod{3^2}\\ \end{aligned} \]

所以:

\[=p_i^{\lfloor \frac {n}{p_i} \rfloor} \times {\lfloor \frac {n}{p_i} \rfloor}! \times (\prod_{i=1,i \not \equiv 0 \ {\rm {mod}} \ p_i} ^{p_i^{k_i}} i)^{\lfloor \frac {n}{p_i^{p_k}} \rfloor} \times \prod_{i=p_i^{k_i} \times \lfloor \frac {n}{p_i^{k_i}} \rfloor,i \not \equiv 0 \ {\rm {mod}} \ p_i} ^{n} i \]

发现其中 \({\lfloor \frac {n}{p_i} \rfloor}!\)\(n!\) 的求解方法是一样的,递归求解即可。

部分参考:Fading大佬的博客


模版:

P4720 【模板】扩展卢卡斯

#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\n",exLucas(n,m,p));
	return 0;
}
posted @ 2020-10-27 15:35  GoPoux  阅读(220)  评论(0编辑  收藏  举报