「笔记」组合数取模
组合数取模
求解：
小 trick
杨辉三角递推。
\(n,m\) 均较小，或 \(n\) 较大且 \(m\) 较小时。
可以得到所有组合数，复杂度都是 \(O(nm)\)。
阶乘 + 逆元。
考虑组合数公式 \(\dfrac{n!}{m!(n-m)!}\)
预处理阶乘，并求逆元即可，复杂度 \(O(n)\) 级别。
以上这两种做法的复杂度与 \(p\) 无关。
Lucas 定理
\(n,m\) 较大，模数 \(p\) 较小且为质数。
描述
有：
发现 \(\left\lfloor\frac{n}{p}\right\rfloor, \left\lfloor\frac{m}{p}\right\rfloor<p\)，这一项可以直接 阶乘+逆元 求得。
要求 \(p\) 不能太大。
递归求解即可，边界为 \(m=0\)，返回 \(1\)。
证明
公式挺好背的，先咕着。
代码
预处理 \(O(p)\)，单次查询 \(O(\log n)\) 级别。
//知识点：Lucas
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstdlib>
#define LL long long
const int kN = 1e5 + 10;
//=============================================================
LL n, m, mod, fac[kN], inv[kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
LL qmul(LL x, LL y) {
LL ret = 1;
while (y) {
if (y & 1) ret = ret * x % mod;
x = x * x % mod, y >>= 1ll;
}
return ret;
}
LL C(int n, int m) {
if (m > n) return 0;
if (m == n) return 1;
return fac[n] * inv[m] % mod * inv[n - m] % mod;
}
LL Lucas(int n, int m) {
return m ? C(n % mod, m % mod) * Lucas(n / mod, m / mod) % mod : 1;
}
//=============================================================
int main() {
fac[0] = 1;
int T = read();
while (T--) {
n = read(), m = read(), mod = read();
for (int i = 1; i <= mod; ++ i) fac[i] = 1ll * fac[i - 1] * i % mod;
inv[mod - 1] = qmul(fac[mod - 1], mod - 2);
for (int i = mod - 2; i >= 1; -- i) inv[i] = 1ll * inv[i + 1] * (i + 1) % mod;
printf("%lld\n", Lucas(n + m, n));
}
return 0;
}
exLucas
和 Lucas 定理并没有关系的一种算法，因用处相同而得名。
并不是一种定理。
\(n,m\) 较大，模数 \(p\) 较小，可不为质数。
算法描述
考虑质因数分解， \(p = \prod\limits_{i=1}p_i^{a_i}\)。
答案即为下列同余方程组的解，中国剩余定理合并即可：
考虑求解单个同余方程的解，有：
由于与 \(p^k\) 不一定互质，不一定能求得 \(m!(n-m)\) 的逆元。
考虑提出 \(p\)，使其互质来求逆元，则有：
\(x,y,z\) 代表质因子 \(p\) 的出现次数。
就可以通过 exgcd 求逆元了。
考虑化简，将上式拆开，先考虑求 \(\dfrac{n!}{p^x}\)。
此式等于 \(n!\) 除去所有 \(p\)，考虑展开 \(n!\) 并提出 \(p\)：
上式在模 \(p\) 意义下展开，则对于 \(\prod\limits_{i=1}^{n}i[p\nmid i]\)，有：
出现了循环节，循环节长度为 \(p^k-1\)。
则有：
现在考虑除去质因子 \(p\)。
仅有 \(p^{\lfloor\frac{n}{p}\rfloor}\) 和 \(\lfloor\frac{n}{p}\rfloor!\) 可能存在质因子 \(p\)，后面那一坨不存在因子 \(p\)。
设 \(\lfloor\frac{n}{p}\rfloor!\) 含质因子 \(p\) 的个数为 \(q\)，则有：
发现这玩意很可递推的样子，设 \(f(n) = \dfrac{n!}{p^x}\)，则有：
就可以递推了。
复杂度约为 \(O(\log_p n)\)，边界 \(f(0)=1\)。
再搞一波 \(p^{x-y-z}\)，发现只需求出指数 \(x-y-z\) 即可。
拆开考虑，先看 \(p^x\)，把上面 \(n!\) 的表达式拿过来：
\(x\) 为 \(n!\) 中质因子 \(p\) 的个数，
显然，将 \(p^{\lfloor\frac{n}{p}\rfloor}\lfloor\frac{n}{p}\rfloor!\) 中质因子 \(p\) 提出，即为 \(p^x\)。
\(p^{\lfloor\frac{n}{p}\rfloor}\) 含有 \({\lfloor\frac{n}{p}\rfloor}\) 个 \(p\)，\(\lfloor\frac{n}{p}\rfloor!\) 不好直接求，还是考虑递推。
设 \(g(n)\) 表示 \(n!\) 中含有质因子 \(p\) 的个数，则有：
复杂度约为 \(O(\log_p n)\)，\(n<p\) 时有边界 \(g(n)=0\)。
最后得到：
合并一下就有答案啦。
复杂度分析
求解单个同余方程，预处理循环节，并递推求 \(f,g\)。
复杂度约为 \(O(p^k + \log^2 n)\)。
合并起来的复杂度约为 \(O(\max\{p\} \log P)\)，\(P\) 是模数，\(p\) 为质因子。
总复杂度大概是 \(O(\sum{p^k + \log^2 n})\) 级别，能过。
感觉比较玄的话可以参考下代码。
代码
//知识点：exLucas
/*
By:Luckyblock
*/
#include <cstdio>
#include <ctype.h>
#include <cstring>
#include <algorithm>
#define ll long long
const int kMaxn = 1010;
//=============================================================
ll a[kMaxn], pk[kMaxn];
//=============================================================
inline ll read() {
ll f = 1, w = 0; char ch = getchar();
for (; !isdigit(ch); ch = getchar()) if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void GetMax(int &fir, int sec) {
if (sec > fir) fir = sec;
}
void GetMin(int &fir, int sec) {
if (sec < fir) fir = sec;
}
ll QuickPow(ll x_, ll y_, ll mod_) {
ll ret = 1;
while (y_) {
if (y_ & 1ll) ret = ret * x_ % mod_;
y_ >>= 1ll, x_ = x_ * x_ % mod_;
}
return ret;
}
ll F(ll n_, ll p_, ll pk_, ll rep_) {
if (! n_) return 1;
ll rep = 1, remain = 1; //循环节 和 剩余部分
rep = QuickPow(rep_, n_ / pk_, pk_);
for (ll i = pk_ * (n_ / pk_); i <= n_; ++ i) {
if (i % p_) remain = remain * (i % pk_) % pk_;
}
return F(n_ / p_, p_, pk_, rep_) * rep % pk_ * remain % pk_;
}
ll G(ll n_, ll p_) {
if (n_ < p_) return 0;
return G(n_ / p_, p_) + (n_ / p_);
}
void exgcd(ll a_, ll b_, ll &x_, ll &y_) {
if (! b_) {
x_ = 1, y_ = 0;
return ;
}
exgcd(b_, a_ % b_, x_, y_);
ll tmp = x_;
x_ = y_, y_ = tmp - a_ / b_ * y_;
}
ll Inv(ll x_, ll p_) {
ll x, y;
exgcd(x_, p_, x, y);
return (x + p_) % p_;
}
ll C(ll n_, ll m_, ll p_, ll pk_) {
ll rep = 1; //预处理循环节
for (ll i = 1; i <= pk_; ++ i) {
if (i % p_) rep = rep * i % pk_;
}
ll fn = F(n_, p_, pk_, rep);
ll fm = Inv(F(m_, p_, pk_, rep), pk_);
ll fnm = Inv(F(n_ - m_, p_, pk_, rep), pk_);
ll powp = QuickPow(p_, G(n_, p_) - G(m_, p_) - G(n_ - m_, p_), pk_);
return fn * fm % pk_ * fnm % pk_ * powp % pk_;
}
ll exLucas(ll n_, ll m_, ll p_) {
ll remain = p_, cnt = 0;
for (ll i = 2; i * i <= p_; ++ i) {
if (remain % i) continue ;
pk[++ cnt] = 1;
while (! (remain % i)) pk[cnt] *= i, remain /= i;
a[cnt] = C(n_, m_, i, pk[cnt]);
}
if (remain != 1) {
pk[++ cnt] = remain;
a[cnt] = C(n_, m_, remain, remain);
}
ll ans = 0;
for (ll i = 1; i <= cnt; ++ i) {
ll m = p_ / pk[i], inv = Inv(m, pk[i]);
ans = (ans + a[i] * m % p_ * inv % p_) % p_;
}
return ans;
}
//=============================================================
int main() {
ll n = read(), m = read(), p = read();
printf("%lld\n", exLucas(n, m, p));
return 0;
}
写在最后
参考资料：
想写群征文活动但是一直没时间啊啊啊啊啊啊
