Jijidawang

组合数取模 合集

目录

0. 介绍

组合数取模，即给定 \(n,m,p\)，求

\[\dbinom nm\bmod p \]

下面按数据范围讨论一些解法：

\(p\) / \(n,m\) \(n,m\le 200\) \(n,m\le 10^6\) \(n,m\le 10^{18}\)
\(p\le 10^6\) 且 \(p\) 为素数 杨辉三角暴力递推 Lucas 定理 Lucas 定理
\(p>10^9\) 且 \(p\) 为素数 杨辉三角暴力递推 预处理阶乘及其逆元 exLucas
\(p=p_1p_2p_3\cdots p_k\) 且 \(p_{1\cdots k}\le 2000\) 为互不相同的素数 杨辉三角暴力递推 Lucas 定理 CRT 合并 Lucas 定理 CRT 合并
\(p\) 无特殊限制 杨辉三角暴力递推 exLucas exLucas

（\(p\) 的第二档本来想表达 \(p\) 远大于 \(n\) 的，但这样 \(n,m\le 10^{18}\) 的档就不可做了 qwq）
（\(p\) 的第三档就是 \(p\) 为 square-free number）
（CRT 指中国剩余定理）

时空复杂度：

做法 预处理时间复杂度 询问时间复杂度 空间复杂度
杨辉三角暴力递推 \(O(n^2)\) \(O(1)\) \(O(n^2)\)
预处理阶乘及其逆元 \(O(n)\) \(O(1)\) \(O(n)\)
Lucas 定理 预处理 \(1\sim p\) 的组合数，取决于所用的算法 \(O(p\log n)\) 预处理 \(1\sim p\) 的组合数，取决于所用的算法
Lucas 定理 CRT 合并 预处理 \(1\sim p\) 的组合数，取决于所用的算法 \(O(P\log n+p\log P)\) 预处理 \(1\sim p\) 的组合数，取决于所用的算法
exLucas 无预处理，\(O(1)\) \(O(\sum p^k\log n(\log n-k)+p\log P)\) \(O(1)\)
多组询问 exLucas（模数相同） \(O(\sum p^k)\) \(O(\log p^k​n)\) \(O(1)\)

在 Lucas 定理 CRT 合并中， \(P\) 是模数，\(p\) 是最大质因子

1. 杨辉三角暴力递推

我们有

\[\dbinom nm=\dbinom{n-1}{m-1}+\dbinom{n-1}m \]

（考虑 \(n\) 里选 \(m\)，讨论选不选第一个即可证明）

\(O(n^2)\) 递推即可 .

2. 预处理阶乘及其逆元

众所周知

\[\dbinom nm=\dfrac{n!}{m!(n-m)!}=n!m!^{-1}(n-m)!^{-1} \]

\(O(n)\) 预处理阶乘及其逆元即可

3. Lucas 定理

Lucas 定理：

形式一

设 \(n,m\) 在 \(p\)（\(p\) 是质数）进制下表示为

\[\overline{n_kn_{k-1}n_{k-2}\dots n_1}\,,\overline{m_km_{k-1}m_{k-2}\dots m_1} \]

其中 \(0\le n_i<p\)，\(0\le m_i<p\) .

则有：

\[\dbinom{n}{m}\equiv \dbinom{n_k}{m_k}\dbinom{n_{k-1}}{m_{k-1}}\cdots\dbinom{n_1}{m_1}\pmod p \]

形式二（简单）

\[C^n_m\equiv C^{n\bmod p}_{m\bmod p}\cdot C^{\lfloor n/p\rfloor}_{\lfloor m/p\rfloor} \]

Lucas 定理的证明

设 \(x\in\mathbb N^+\)，且 \(x<p\)，则：

\[C^x_p=\dfrac{p!}{x!(p-x)!}=\dfrac{p\cdot (p-1)!}{x\cdot(x-1)!\cdot(p-x)!}=\dfrac px\cdot C^{x-1}_{p-1} \]

由于 \(x<p\) 且 \(p\) 是素数，所以 \(x\) 存在模 \(p\) 意义下的逆元，因此：

\[C^x_p\equiv p\cdot x^{-1}\cdot C^{x-1}_{p-1}\pmod p \]

因为 \(p\cdot x^{-1}\cdot C^{x-1}_{p-1}\equiv 0\pmod p\) ，故 \(C^x_p\equiv0\pmod p\).

由二项式定理，

\[(1+x)^p=\sum_{i=0}^pC_p^ix_i \]

除了 \(i=0,i=p\) 的项数，别的模 \(p\) 均为 \(0\)，所以：

\[(1+x)^p\equiv 1+x^p\pmod p \]

现在我们设：

\[\begin{cases}\lfloor m/p\rfloor=q_m\\\lfloor n/p\rfloor=q_n\end{cases}\ ,\ \begin{cases}m\bmod p=r_m\\n\bmod p=r_n\end{cases} \]

从而：

\[\begin{cases}m=q_mp+r_m\\n=q_np+r_n\end{cases} \]

再由二项式定理有

\[(1+x)^m=\sum_{i=0}^mC_m^ix_i \]

而同时又有：

\[\begin{aligned}(1+x)^m&=(1+x)^{q_mp+r_m}\\&=(1+x)^{q_mp}\cdot(q+x)^{r_m}\\&=((1+x)^q)^{m_p}\cdot(q+x)^{r_m}\\&\equiv (1+x^p)^{q_m}\cdot(1+x)^{r_m}\\&=\sum_{i=0}^{q_m}C^i_{q_m}x^{ip}\sum_{i=0}^{r_m}C^i_{r_m}x^{i}\\&=\sum_{i=0}^{q_m}\sum_{j=0}^{r_m}C^i_{q_m}C^j_{r_m}x^{ip+j}\pmod p\end{aligned} \]

注意到 \(ip+j\) 正好能取到 \(0\) 到 \(m\) 内所有整数一次，所以枚举 \(k=ip+j\)，得

\[(x+1)^m\equiv \sum_{k=0}^n C_{q_m}^{\lfloor k/p\rfloor}C_{r_m}^{k\bmod p}x^k\pmod p \]

结合二项式定理，得

\[\sum_{k=0}^mC^k_mx^k\equiv \sum_{k=0}^n C_{q_m}^{\lfloor k/p\rfloor}C_{r_m}^{k\bmod p}x^k\pmod p \]

对比系数，得：

\[C^n_m\equiv C^{n\bmod p}_{m\bmod p}\cdot C^{\lfloor n/p\rfloor}_{\lfloor m/p\rfloor}=C_{q_m}^{q_n}C_{r_m}^{r_n}\pmod p \]

命题获证 .

用形式一可以写出一个递推的代码：

int Lucas(int n,int m,int p)
{
	for (int i=0;i<p;i++)
	{
		C[i][0]=1;
		for (int j=1;j<=i;j++)
		{
			C[i][j]=C[i-1][j-1]+C[i-1][j];
			if (C[i][j]>=p) C[i][j]-=p;
		}
	} // init
	int ans=1;
	while (n||m){ans=ans*C[n%p][m%p]; n/=p; m/=p;}
	return ans;
}

而用形式二则可以写出一个递归的代码：

ll lucas(ll n,ll m){return (n<m)?0:(n<p)?C(n,m):C(n%p,m%p)*lucas(n/p,m/p)%p;}

4. Lucas 定理 CRT 合并

令

\[A=\dbinom nm \]

则

\[\begin{cases}A\equiv \dbinom nm\pmod{p_1}\\A\equiv \dbinom nm\pmod{p_2}\\\cdots\\A\equiv \dbinom nm\pmod{p_k}\end{cases} \]

用 Lucas 定理求组合数，然后 CRT 合并即可

正确性显然

5. exLucas

先把 \(p\) 分解质因数，令

\[p=\prod_{i=1}^t p_i^{k_i} \]

用同样的做法，列出同余方程组，求出 \(\dbinom nm\bmod p_i^{k_i}\) 然后 CRT 合并 .

那么问题来了，\(\dbinom nm\bmod p_i^{k_i}\) 怎么求？

我们知道，\(\dbinom nm=\dfrac{n!}{m!(n-m)!}\)，由于逆元可能不存在，考虑提取出质因子 \(p_i\) 来

\[\dfrac{n!}{m!(n-m)!}\equiv \dfrac{\dfrac{n!}{p_i^A}}{\dfrac{m!}{p_i^B}\cdot \dfrac{(n-m)!}{p_i^C}}\cdot p_i^{A-B-C}\pmod{p_i^{k_i}} \]

其中 \(A\) 为 \(n!\) 中质因子 \(p_i\) 的次数，\(B,C\) 同理 .

此时分母就存在逆元了，现在问题又转化成了求形如

\[\dfrac{n!}{p^a}\bmod p_k \]

先考虑如何计算 \(n!\bmod p_k\) .

举个例子（\(22!\bmod 3^2\)）

\[\begin{aligned}22!&\equiv1\times2\times3\times4\times5\times6\times7\times8\times9\times10\times11\times12\times13\times14\times15\times16\times17\times18\times19\times20\times21\times22\\&\equiv3^7\times(1\times2\times3\times4\times5\times6\times7)\times(1\times2\times4\times5\times7\times8\times10\times11\times13\times14\times16\times17\times19\times20\times22)\pmod{p_k}\end{aligned} \]

可以看出上式分为三个部分：
第一部分是 \(p\) 的幂，次数是小于等于 \(n\) 的 \(p\) 的倍数的个数（即 \(\left\lfloor\dfrac{n}{p}\right\rfloor\)）.

第二部分是一个阶乘 \(7!\)，即 \(\left\lfloor\dfrac{n}{p}\right\rfloor!\)，可以递归求

第三部分是 \(n!\) 中与 \(p\) 互质的部分的乘积，这一部分具有如下性质（加个 \(p^k\)）：

\[1\times 2\times 4\times 5\times 7\times 8\equiv10\times 11\times 13\times 14\times 16\times 17\pmod{p^k} \]

一般的，有

\[\forall i,\prod_{i,(i,p)=1}^{p^k} i\equiv \prod_{i,(i,p)=1}^{p^k}(i+tp^k)\pmod{p^k} \]

整块快速幂小块暴力即可求出 \(n!\bmod p^k\)

回到原问题，求

\[\dfrac{n!}{p^a}\bmod p_k \]

分母全部由上述第一部分和第二部分贡献（第三部分和 \(p\) 互质）。而递归计算第二部分的时候已经除去了第二部分中的因子 \(p\)，所以最终的答案就是上述第二部分递归返回的结果和第三部分的乘积（与第一部分无关） .

然后一层层回去即可求出组合数取模 .

typedef long long ll;
ll qpow(ll a,ll n,ll p)
{
	ll ans=1;
	while (n)
	{
		if (n&1) ans=ans*a%p;
		a=a*a%p; n>>=1;
	} return ans;
}
namespace extra_lucas_detail
{
	ll fac(ll n,ll p,ll pk)
	{
		if (!n) return 1;
		ll ans=1;
		for (int i=1;i<pk;i++)
			if (i%p) ans=ans*i%pk;
		ans=qpow(ans,n/pk,pk);
		for (int i=1; i<=n%pk;i++)
			if (i%p) ans=ans*i%pk;
		return ans*fac(n/p,p,pk)%pk;
	}
	ll exgcd(ll a,ll b,ll& x,ll& y)
	{
	    if (!b){x=1; y=0; return a;}
	    ll d=exgcd(b,a%b,x,y);
		int z=x; x=y; y=z-y*(a/b);
		return d;
	}
	ll inv(ll a,ll P){ll x,y,d=exgcd(a,P,x,y); return (d==1)?x%P:-1;}
	ll C(ll n,ll m,ll p,ll pk)
	{
		if (n<m) return 0;
		ll A=fac(n,p,pk),B=fac(m,p,pk),C=fac(n-m,p,pk),cnt=0;
		for (ll i=n;i;i/=p) cnt+=i/p;
		for (ll i=m;i;i/=p) cnt-=i/p;
		for (ll i=n-m;i;i/=p) cnt-=i/p;
		return A*inv(B,pk)%pk*inv(C,pk)%pk*qpow(p,cnt,pk)%pk;
	}
}
ll exlucas(ll n,ll m,ll p)
{
	using namespace extra_lucas_detail;
	ll ans=0,P=p;
	for (int i=2;(p>1)&&(i*i<=p);i++)
	{
		ll pk=1;
		while (!(p%i)){p/=i; pk*=i;}
		if (pk>1)
		{
			ll now=C(n,m,i,pk);
			ans=(ans+now*(P/pk)%P*inv(P/pk,pk)%P)%P;
		}
	}
	if (p>1)
	{
		ll now=C(n,m,p,p);
		ans=(ans+now*(P/p)%P*inv(P/p,p)%P)%P;
	} return (ans%P+P)%P;
}
posted @ 2021-07-01 17:18  yspm  阅读(228)  评论(0)    收藏  举报
😅​