Description

求

\[\dbinom{n}{m}\bmod p \]

其中 \(p\) 較小且 不保證 \(p\) 是質數。

Method

前置芝士：

• 中國剩餘定理

因為 \(p\) 不為質數，所以得使用擴充套件盧卡斯定理（\(\rm Extended\ Lucas,exLucas\)）。

Part 1：中國剩餘定理

首先按照 唯一分解定理 將 \(p\) 分解質因數：

\[p=\prod\limits_{i=1}^{k}q_i^{\alpha_i} \]

如果我們能夠求出

\[\begin{cases} \dbinom{n}{m}\bmod q_1^{\alpha_1}\\ \dbinom{n}{m}\bmod q_2^{\alpha_2}\\ \cdots\\ \dbinom{n}{m}\bmod q_k^{\alpha_k} \end{cases} \]

那麼因為 \(\forall i\ne j,\gcd(q_i^{\alpha_i},q_j^{\alpha_j})=1\)

，所以可以用 CRT 合併。

那麼問題就變成了求出 \(\dbinom{n}{m}\bmod q_i^{\alpha_i}\)。

Part 2：移除分子分母中的質因數

\[\begin{aligned} \dbinom{n}{m} &\equiv \dfrac{n!}{m!(n-m)!}\\ &\equiv n!\cdot \operatorname{inv}(m!)\cdot \operatorname{inv}[(n-m)!] \pmod {q^{\alpha}} \end{aligned} \]

但 \(m!,(n-m)!\) 與 \(q^{\alpha}\) 不一定互質，所以不一定存在逆元。

我們將分子和分母中的 \(q\) 的倍數全部提出來，設 \(n!,m!,(n-m)!\) 分別含 \(x,y,z\) 個質因數 \(q\)。

\[\begin{aligned} \dfrac{n!}{m!(n-m)!} &=\dfrac{\dfrac{n!}{q^x}}{\dfrac{m!}{q^y}\cdot \dfrac{(n-m)!}{q^z}}\cdot q^{x-y-z}\\ &\equiv \dfrac{n!}{q^x}\cdot \operatorname{inv}\left(\dfrac{m!}{q^y}\right)\cdot \operatorname{inv}\left[\dfrac{(n-m)!}{q^z}\right]\pmod {q^{\alpha}} \end{aligned} \]

這時 \(\dfrac{m!}{q^y},\dfrac{(n-m)!}{q^z}\)

都與 \(q^{\alpha}\) 互質。

Part 3：計算

以 \(n!\) 為例：

\[\begin{aligned} n! &=(q\times 2q\times \cdots\times \left\lfloor\dfrac{n}{q}\right\rfloor q)\times \prod\limits_{q\nmid i}^n i\\ &=q^{\left\lfloor\frac{n}{q}\right\rfloor}\times \left\lfloor\dfrac{n}{q}\right\rfloor!\times \prod\limits_{q\nmid i}^n i \end{aligned} \]

其中 \(q^{\left\lfloor\frac{n}{q}\right\rfloor}\) 可用快速冪直接算，\(\left\lfloor\dfrac{n}{q}\right\rfloor!\) 可遞迴求解，後面 \(\prod\limits_{q\nmid i}^n i\) 就不好算了。

我們發現因為模數是 \(q^{\alpha}\)，又有 \(x+q^{\alpha}\equiv x\pmod{q^{\alpha}}\)，所以可以將 \(q^{\alpha}\) 作為一個迴圈週期，暴力算出 \(\prod\limits_{q\nmid i}^{q^{\alpha}}i\)，然後迴圈實際上有 \(\left\lfloor\dfrac{n}{q^{\alpha}}\right\rfloor\) 個，直接用快速冪。

對於多出來的部分，長度一定小於 \(q^{\alpha}\)，也可以暴力算。

以 \(n=19,p=3,\alpha=2\) 為例：

\[\begin{aligned} 19! &=1\times2\times3\times4\times5\times6\times7\times8\times9\times10\times11\times12\times13\times14\times15\times16\times17\times18\times19\\ &=(1\times2\times4\times5\times7\times8)\times(10\times11\times13\times14\times16\times17)\times19\times(3\times6\times9\times12\times15\times18)\\ &=(1\times2\times4\times5\times7\times8)\times(10\times11\times13\times14\times16\times17)\times19\times3\times(1\times2\times3\times4\times5\times6)\\ &\equiv (1\times2\times4\times5\times7\times8)^2\times19\times3^6\times6!\pmod {3^2}\\ &=\prod\limits_{3\nmid i}^{19}i\times 3^{\left\lfloor\frac{19}{3}\right\rfloor}\times \left\lfloor\dfrac{19}{3}\right\rfloor!~ \end{aligned} \]

再遞迴計算 \(6!\bmod 3^2\) 即可。

遞迴部分的程式碼如下：

ll cal(ll n, ll p, ll pa)
{
	if (!n)
	{
		return 1;
	}
	ll ans = 1;
	for (int i = 1; i <= pa; i++) //迴圈部分
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	ans = qpow(ans, n / pa, pa);
	for (int i = 1; i <= n % pa; i++) //多出部分
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	return ans * cal(n / p, p, pa) % pa;
}

注意到這部分其實求的就是 \(\dfrac{n!}{q^x}\)。

對於 \(m!,(n-m)!\) 的處理同理。

Part 4：合併

通過 Part 3 可以求出 \(\dfrac{n!}{q^x},\dfrac{m!}{q^y},\dfrac{(n-m)!}{q^z}\)，現在只需要求出 \(x,y,z\) 就可以得到 \(q^{x-y-z}\) 了。

怎麼算想必小奧都講過了吧，比如說

\[x=\left\lfloor\dfrac{n}{q}\right\rfloor+\left\lfloor\dfrac{n}{q^2}\right\rfloor+\left\lfloor\dfrac{n}{q^3}\right\rfloor+\cdots \]

至此可求出 \(\dbinom{n}{m}\bmod q_i^{\alpha_i}\)，再結合 Part 1 使用 CRT 合併就完成了求解 \(\dbinom{n}{m}\bmod p\) 的過程。

Code

//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

ll qpow(ll a, ll b, ll p)
{
	ll base = a, ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = ans * base % p;
		}
		base = base * base % p;
		b >>= 1;
	}
	return ans;
}

ll cal(ll n, ll p, ll pa)
{
	if (!n)
	{
		return 1;
	}
	ll ans = 1;
	for (int i = 1; i <= pa; i++)
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	ans = qpow(ans, n / pa, pa);
	for (int i = 1; i <= n % pa; i++)
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	return ans * cal(n / p, p, pa) % pa;
}

ll cnt_p(ll n, ll m, ll p)
{
	ll cnt = 0;
	for (ll i = p; i <= n; i *= p)
	{
		cnt += n / i;
	}
	for (ll i = p; i <= m; i *= p)
	{
		cnt -= m / i;
	}
	for (ll i = p; i <= n - m; i *= p)
	{
		cnt -= (n - m) / i;
	}
	return cnt;
}

ll x, y;

void exgcd(ll a, ll b)
{
	if (!b)
	{
		x = 1, y = 0;
		return;
	}
	exgcd(b, a % b);
	ll tmp = x;
	x = y;
	y = tmp - a / b * y;
}

ll inv(ll a, ll p)
{
	exgcd(a, p);
	x = (x % p + p) % p;
	return x;
}

ll C(ll n, ll m, ll p, ll pa)
{
	ll a = cal(n, p, pa), b = cal(m, p, pa), c = cal(n - m, p, pa), cnt = cnt_p(n, m, p);
	return a * inv(b, pa) % pa * inv(c, pa) % pa * qpow(p, cnt, pa) % pa;
}

ll a[10], b[10];

ll CRT(int n)
{
	ll m = 1;
	for (int i = 1; i <= n; i++)
	{
		m *= a[i];
	}
	ll ans = 0;
	for (int i = 1; i <= n; i++)
	{
		ll mi = m / a[i];
		ll Mi = inv(mi, a[i]);
		ans = (ans + b[i] * mi % m * Mi % m) % m;
	}
	return ans;
}

ll exLucas(ll n, ll m, ll p)
{
	int k = 0;
	for (ll i = 2; i * i <= p; i++)
	{
		if (p % i == 0)
		{
			a[++k] = 1;
			while (p % i == 0)
			{
				a[k] *= i;
				p /= i;
			}
			b[k] = C(n, m, i, a[k]);
		}
	}
	if (p > 1)
	{
		a[++k] = p;
		b[k] = C(n, m, p, p);
	}
	return CRT(k);
}

int main()
{
	ll n, m, p;
	scanf("%lld%lld%lld", &n, &m, &p);
	printf("%lld\n", exLucas(n, m, p));
	return 0;
}

Complexity

計算 \(\dfrac{n!}{q^x}\bmod q^{\alpha}\) 的時間為 \(O(q^{\alpha}\log n)\)。

計算 \(x\) 的時間為 \(O(\log n)\)。

所以計算 \(\dbinom{n}{m}\bmod q^{\alpha}\) 的時間為 \(O(q^{\alpha}\log n)\)。

整個預處理部分就是 \(O(\sqrt{p}+\sum\limits_{i=1}^k q_i^{\alpha_i}\log n)\)。

CRT 部分是 \(O(k\log \prod\limits_{i=1}^k q_i^{\alpha_i})=O(k\log p)\) 的。

綜上，exLucas 的時間複雜度為 \(O(\sqrt{p}+\sum\limits_{i=1}^k q_i^{\alpha_i} \log n+k\log p)\sim O(p\log p)\)。

Problems

模板

P4720 【模板】擴充套件盧卡斯定理/exLucas

模板 + 乘法原理

P2183 [國家集訓隊]禮物

答案為 \(\dbinom{n}{w_1}\times \dbinom{n-w_1}{w_2}\times \dbinom{n-w_1-w_2}{w_3}\times \cdots\bmod P\)。

References

• [1] JokerJim：【演算法】擴充套件盧卡斯詳解
• [2] YangTY：Lucas 定理
• [3] OI Wiki：盧卡斯定理