「笔记」组合数取模


目录
  • 组合数取模
  • 小 trick
  • Lucas 定理
    • 描述
    • 证明
    • 代码
  • exLucas
    • 算法描述
  • 写在最后


组合数取模

求解:

\[{n\choose m} \mod p \]


小 trick

杨辉三角递推。
\(n,m\) 均较小,或 \(n\) 较大且 \(m\) 较小时。
可以得到所有组合数,复杂度都是 \(O(nm)\)

阶乘 + 逆元。
考虑组合数公式 \(\dfrac{n!}{m!(n-m)!}\)
预处理阶乘,并求逆元即可,复杂度 \(O(n)\) 级别。

以上这两种做法的复杂度与 \(p\) 无关。


Lucas 定理

\(n,m\) 较大,模数 \(p\) 较小且为质数。

描述

有:

\[{n\choose m} \bmod p = {\lfloor\frac{n}{p}\rfloor\choose \lfloor\frac{m}{p}\rfloor}\cdot {n\bmod p\choose m\bmod p} \bmod p \]

发现 \(\left\lfloor\frac{n}{p}\right\rfloor, \left\lfloor\frac{m}{p}\right\rfloor,这一项可以直接 阶乘+逆元 求得。
要求 \(p\) 不能太大。
递归求解即可,边界为 \(m=0\),返回 \(1\)

证明

公式挺好背的,先咕着。

代码

P3807 【模板】卢卡斯定理

预处理 \(O(p)\),单次查询 \(O(\log n)\) 级别。

//知识点:Lucas
/*
By:Luckyblock
*/
#include 
#include 
#include 
#include 
#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}\)
答案即为下列同余方程组的解,中国剩余定理合并即可:

\[\begin{cases} {n\choose m} \pmod{p_1^{a_1}}\\ {n\choose m} \pmod{p_2^{a_2}}\\ \cdots\\ {n\choose m} \pmod{p_i^{a_i}} \end{cases}\]


考虑求解单个同余方程的解,有:

\[\begin{aligned} {n\choose m}\equiv\dfrac{n!}{m!(n-m)!}\pmod {p^k} \end{aligned}\]

由于与 \(p^k\) 不一定互质,不一定能求得 \(m!(n-m)\) 的逆元。
考虑提出 \(p\),使其互质来求逆元,则有:

\[\begin{aligned} &\dfrac{n!}{m!(n-m)!}\equiv \dfrac{\dfrac{n!}{p^x}}{\dfrac{m!}{p^y}\dfrac{(n-m)!}{p^z}}p^{x-y-z}\pmod {p^k} \end{aligned}\]

\(x,y,z\) 代表质因子 \(p\) 的出现次数。
就可以通过 exgcd 求逆元了。


考虑化简,将上式拆开,先考虑求 \(\dfrac{n!}{p^x}\)
此式等于 \(n!\) 除去所有 \(p\),考虑展开 \(n!\) 并提出 \(p\)

\[\begin{aligned} n! &= 1\times 2 \times 3\cdots \times n\\ &= (p \times 2p\times \cdots\times \left\lfloor\dfrac{n}{p}\right\rfloor p)(1\times 2\times \cdots)\\ &= p^{\lfloor\frac{n}{p}\rfloor}\lfloor\frac{n}{p}\rfloor! \prod_{i=1}^{n}i[p\nmid i] \end{aligned}\]

上式在模 \(p\) 意义下展开,则对于 \(\prod\limits_{i=1}^{n}i[p\nmid i]\),有:

\[\begin{aligned} \prod\limits_{i=1}^{n}i[p\nmid i] &= \prod\limits_{i=1}^{n}(i\bmod p^k)[p\nmid i]\\ &= \left(\prod\limits_{i=1}^{p^k-1}i[p\nmid i]\right)^{\left\lfloor\frac{n}{p^k}\right\rfloor}\left(\prod_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor}^ni[p\nmid i]\right) \end{aligned}\]

出现了循环节,循环节长度为 \(p^k-1\)
则有:

\[n! = p^{\lfloor\frac{n}{p}\rfloor}\lfloor\frac{n}{p}\rfloor! \left(\prod\limits_{i=1}^{p^k}i[p\nmid i]\right)^{\left\lfloor\frac{n}{p^k}\right\rfloor}\left(\prod_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor}^ni[p\nmid i]\right) \]

现在考虑除去质因子 \(p\)
仅有 \(p^{\lfloor\frac{n}{p}\rfloor}\)\(\lfloor\frac{n}{p}\rfloor!\) 可能存在质因子 \(p\),后面那一坨不存在因子 \(p\)
\(\lfloor\frac{n}{p}\rfloor!\) 含质因子 \(p\) 的个数为 \(q\),则有:

\[\dfrac{n!}{p^x} = \dfrac{\lfloor\frac{n}{p}\rfloor!}{p^q}\left(\prod\limits_{i=1}^{p^k}i[p\nmid i]\right)^{\left\lfloor\frac{n}{p^k}\right\rfloor}\left(\prod_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor}^ni[p\nmid i]\right) \]

发现这玩意很可递推的样子,设 \(f(n) = \dfrac{n!}{p^x}\),则有:

\[f(n) = f(\left\lfloor\frac{n}{p}\right\rfloor)\left(\prod\limits_{i=1}^{p^k}i[p\nmid i]\right)^{\left\lfloor\frac{n}{p^k}\right\rfloor}\left(\prod_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor}^ni[p\nmid i]\right) \]

就可以递推了。
复杂度约为 \(O(\log_p n)\),边界 \(f(0)=1\)


再搞一波 \(p^{x-y-z}\),发现只需求出指数 \(x-y-z\) 即可。

拆开考虑,先看 \(p^x\),把上面 \(n!\) 的表达式拿过来:

\[n! = p^{\lfloor\frac{n}{p}\rfloor}\lfloor\frac{n}{p}\rfloor! \left(\prod\limits_{i=1}^{p^k}i[p\nmid i]\right)^{\left\lfloor\frac{n}{p^k}\right\rfloor}\left(\prod_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor}^ni[p\nmid i]\right) \]

\(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\) 的个数,则有:

\[g(n)=\lfloor\frac{n}{p}\rfloor + g(\left\lfloor\frac{n}{p}\right\rfloor) \]

复杂度约为 \(O(\log_p n)\)\(n 时有边界 \(g(n)=0\)


最后得到:

\[{n\choose m} \equiv \dfrac{f(n)}{f(m)f(n-m)}p^{g(n)-g(m)-g(n-m)} \pmod {p^k} \]

合并一下就有答案啦。


复杂度分析

求解单个同余方程,预处理循环节,并递推求 \(f,g\)
复杂度约为 \(O(p^k + \log^2 n)\)
合并起来的复杂度约为 \(O(\max\{p\} \log P)\)\(P\) 是模数,\(p\) 为质因子。

总复杂度大概是 \(O(\sum{p^k + \log^2 n})\) 级别,能过。
感觉比较玄的话可以参考下代码。


代码

P4720 【模板】扩展卢卡斯

//知识点:exLucas
/*
By:Luckyblock
*/
#include 
#include 
#include 
#include 
#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;
}

写在最后

参考资料:

P4720 【模板】扩展卢卡斯 Fading 的题解

想写群征文活动但是一直没时间啊啊啊啊啊啊