【数学】扩展卢卡斯定理
Description
求
其中 \(p\) 较小且 不保证 \(p\) 是质数。
Method
前置芝士：
- 中国剩余定理
因为 \(p\) 不为质数，所以得使用扩展卢卡斯定理（\(\rm Extended\ Lucas,exLucas\)）。
Part 1：中国剩余定理
首先按照 唯一分解定理 将 \(p\) 分解质因数：
如果我们能够求出
那么因为 \(\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：移除分子分母中的质因数
但 \(m!,(n-m)!\) 与 \(q^{\alpha}\) 不一定互质，所以不一定存在逆元。
我们将分子和分母中的 \(q\) 的倍数全部提出来，设 \(n!,m!,(n-m)!\) 分别含 \(x,y,z\) 个质因数 \(q\)。
这时 \(\dfrac{m!}{q^y},\dfrac{(n-m)!}{q^z}\) 都与 \(q^{\alpha}\) 互质。
Part 3：计算
以 \(n!\) 为例：
其中 \(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\) 为例：
再递归计算 \(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}\) 了。
怎么算想必小奥都讲过了吧，比如说
至此可求出 \(\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 n)\)。
Problems
模板
模板 + 乘法原理
答案为 \(\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：卢卡斯定理
