(扩展)Lucas定理——求大组合数取模

最新推荐文章于 2022-02-04 18:33:29 发布

CaptainHarryChen

最新推荐文章于 2022-02-04 18:33:29 发布

阅读量336

点赞数 1

CC 4.0 BY-SA版权

分类专栏： OI知识总结计算机数学

本文链接：https://blog.youkuaiyun.com/can919/article/details/79765790

OI知识总结同时被 2 个专栏收录

37 篇文章

订阅专栏

计算机数学

29 篇文章

订阅专栏

本文介绍了Lucas定理及其扩展版本，详细解释了如何利用这些定理解决组合数取模的问题，包括质数和非质数情况。通过具体实例展示了如何计算阶乘并提取质因子，最终求解复杂场景下的组合数。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Lucas定理

若 $P$ 为质数
$P)C_n^r\equiv C_{\lfloor n/P\rfloor}^{\lfloor r/P\rfloor}\times C_{n\space mod\space P}^{r\space mod\space P}\space (mod \space P)$
要用生成函数证明（不会，所以不证了）

实现直接递归：

void exgcd(long long a,long long b,long long &x,long long &y)
{
	if(b==0)
		x=1,y=0;
	else
	{
		long long x1,y1;
		exgcd(b,a%b,x1,y1);
		x=y1;
		y=x1-(a/b)*y1;
	}
}
long long inv(long long x,long long p)
{
	long long i,t;
	exgcd(x,p,i,t);
	i=(i%p+p)%p;
	return i;
}


long long C(long long n,long long r,long long MOD)
{
	if(n<r)
		return 0;
	if(r==0)
		return 1;
	long long ans=1LL;
	ans=C(n/MOD,r/MOD,MOD);
	n%=MOD,r%=MOD;
	long long a=1,b=1;
	for(int i=n;i>n-r;i--)
		a=(a*i)%MOD;
	for(int i=1;i<=r;i++)
		b=(b*i)%MOD;
	ans=(ans*a%MOD*inv(b,MOD))%MOD;
	return ans;
}

扩展Lucas定理

可以解决模数 $P$ 不是质数的组合数取模

将 $P$ 质因数分解： $P=p_1^{k_1}p_2^{k_2}...p_t^{k_t}$
如果我们求出一下方程中的 $a_1,a_2...a_t$ ，就可以利用中国剩余定理求出 $C_n^r$
$\left \{ \begin{array}{c} C_n^r\equiv a_1\space (mod\space p_1^{k_1})\\ C_n^r\equiv a_2\space (mod\space p_2^{k_2})\\ C_n^r\equiv a_3\space (mod\space p_3^{k_3})\\ ...\\ C_n^r\equiv a_t\space (mod\space p_t^{k_t})\\ \end{array} \right .$

现在来求 $pk)C_n^r\equiv a\space(mod\space p^k)$
$pk)C_n^r\equiv \frac {n!} {(n-r)!\space r!}\space (mod\space p^k)$
令 $pk)a\equiv n!\space(mod\space p^k)$ ， $pk)b\equiv r!\space(mod\space p^k)$ ， $pk)c\equiv(n-r)!\space(mod\space p^k)$

将组合数分成三个阶乘计算： $pk)C_n^r\equiv a\space b^{-1}c^{-1}\space(mod\space p^k)$
（ $b^{-1}$ 表示逆元）
但此时如果 $b$ 与 $p^k$ 不互质，就无法求出逆元，所以需要把 $a$ , $b$ , $c$ 中的质因子 $p$ 全部提出来，求了逆元之后再乘回去。

举个例子：计算 $\space(mod\space 3^2)$
$\begin{aligned} 19!&=1\times 2\times 3\times 4\times ...\times 19\\ &= (1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19)\times3^6 \times (1\times 2\times 3\times 4\times 5\times 6)\\ \end{aligned}$
即把阶乘中的 $p$ 倍数提出来，后面 $6!$ 继续递归求解；
第一个括号里的式子有周期性的：
$32)1\times 2\times 4\times 5\times 7\times 8\equiv 10\times 11\times 13\times 14\times 16\times 17 \space (mod \space 3^2)$
所以暴力求出一个区间的值 $t$ ，结果乘 $t⌊199⌋t^{\lfloor \frac {19} 9 \rfloor}$

细节见代码：

long long pow_mod(long long a,long long b,long long p)
{
	long long res=1;
	while(b)
	{
		if(b&1LL)
			res=(res*a)%p;
		b>>=1LL;
		a=(a*a)%p;
	}
	return res;
}
void exgcd(long long a,long long b,long long &x,long long &y)
{
	if(b==0)
	{
		x=1,y=0;
		return;
	}
	exgcd(b,a%b,y,x);
	y-=(a/b)*x;
}
long long inv(long long a,long long p)//逆元必须用扩展欧几里得，因为p不为质数
{
	long long i,t;
	exgcd(a,p,i,t);
	i=(i%p+p)%p;
	return i;
}

long long Fac(long long N,long long p,long long pk)
{
	if(N==0)
		return 1;
	long long ans=1;
	if(N>=pk)//周期性
	{
		for(long long i=1;i<pk;i++)
			if(i%p)
				ans=(ans*i)%pk;
		ans=pow_mod(ans,N/pk,pk);
	}
	for(long long i=1;i<=N%pk;i++)//除了周期性的，剩下的几个数
		if(i%p)
			ans=(ans*i)%pk;
	ans=(ans*Fac(N/p,p,pk))%pk;//递归求解
	return ans;
}
long long C(long long N,long long R,long long p,long long pk)
{
	if(R==0)
		return 1;
	long long a=Fac(N,p,pk),b=Fac(N-R,p,pk),c=Fac(R,p,pk),k=0,ans;//求出几个阶乘，不包含质因子p
	//下面单独计算质因子p的个数
	for(long long i=N;i;i/=p)
		k+=i/p;
	for(long long i=N-R;i;i/=p)
		k-=i/p;
	for(long long i=R;i;i/=p)
		k-=i/p;
	ans=a*inv(b,pk)%pk*inv(c,pk)%pk*pow_mod(p,k,pk)%pk;
	return ans;
}
long long C(long long N,long long R,long long M)
{
	if(N<R)
		return 0;
	long long ans=0,m=M,temp;
	for(long long p=2,pk;p*p<=m;p++)
		if(m%p==0)//枚举模数M的因子p
		{
			pk=1;
			while(m%p==0)
				m/=p,pk*=p;
			temp=C(N,R,p,pk);//求C(N,R)mod p^k
			ans=(ans+temp*(M/pk)%M*inv(M/pk,pk)%M)%M;//中国剩余定理
		}
	if(m!=1)
	{
		temp=C(N,R,m,m);
		ans=(ans+temp*(M/m)%M*inv(M/m,m)%M)%M;
	}
	return ans;
}