扩展卢卡斯定理模板 / p4720

最新推荐文章于 2024-05-20 20:08:58 发布

Lusfiee

最新推荐文章于 2024-05-20 20:08:58 发布

阅读量175

点赞数

分类专栏：数学习题文章标签：算法

本文链接：https://blog.csdn.net/m0_54646258/article/details/125810186

版权

数学习题专栏收录该内容

13 篇文章 0 订阅

订阅专栏

文章目录

前言
一、题目
二、思路及代码
- 1.思路
- 2.代码

前言

扩展卢卡斯算法可以解决 $C_n^m \% mod$ ，mod 不是素数的问题

可以在 $O(mod(\log mod))$ 时间内处理 $1 e 18$ 的 m 和 n

一、题目

题目链接洛谷P4720

二、思路及代码

1.思路

可以利用中国剩余定理，只需求 $\dbinom{m}{n}\equiv r_i \pmod{p^k}$ 即可
利用威尔逊定理和逆元、快速幂的知识，转化为求：
$\displaystyle \frac{\frac{n!}{q^x}}{\frac{m!}{q^y}\frac{(n-m)!}{q^z}}q^{x-y-z} \pmod{q^k}$
不是很好理解，留下板子即可

2.代码

代码如下：

#include <iostream>
#define int long long
using namespace std;
int n, m, mod;
int quickpow(int base, int x, int p) {
  int ans = 1;
  while (x) {
    if (x & 1) ans = ans * base % p;
    base = base * base % p;
    x >>= 1;
  }
  return ans;
}
int exgcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  }
  int r = exgcd(b, a % b, x, y);
  int t = x;
  x = y;
  y = t - a / b * y;
  return r;
}
int inv(int a, int b) {
  int x, y;
  return exgcd(a, b, x, y) == 1 ? (x + b) % b : -1;
}
int CRT(int x, int ps) { return x * inv(mod / ps, ps) % mod * mod / ps % mod; }
int fac(int x, int ps, int p) {  // ps = p ^ k
  if (!x) return 1;
  int res = 1;
  for (int i = 2; i <= ps; i++)
    if (i % p) res = res * i % ps;
  res = quickpow(res, x / ps, ps);
  for (int i = 2; i <= x % ps; i++)
    if (i % p) res = res * i % ps;
  return res * fac(x / p, ps, p) % ps;  // 不含 p 的 x! % ps
}
int Comp(int n, int m, int ps, int p) {  // ps = p ^ k, 返回组合数模 ps 的值
  int fn = fac(n, ps, p), fm = fac(m, ps, p), fnm = fac(n - m, ps, p);
  int k = 0;
  for (int i = n; i; i /= p) k += i / p;
  for (int i = m; i; i /= p) k -= i / p;
  for (int i = n - m; i; i /= p) k -= i / p;  // 推组合数中 p 的幂次(勒让德)
  return fn * inv(fm, ps) % ps * inv(fnm, ps) % ps * quickpow(p, k, ps) % ps;
}
int exlucas(int n, int m) {  // C_n^m % mod
  int q = mod, ans = 0;
  for (int i = 2; i * i <= mod; i++)
    if (q % i == 0) {
      int ps = 1;
      while (q % i == 0) ps *= i, q /= i;  // ps * q = mod, ps = i(p) ^ k
      ans = (ans + CRT(Comp(n, m, ps, i), ps)) % mod;
    }
  if (q > 1) ans = (ans + CRT(Comp(n, m, q, q), q) % mod);
  return ans % mod;
}
signed main() {
  //   freopen("in.txt", "r", stdin);
  //   freopen("out.txt", "w", stdout);
  scanf("%lld%lld%lld", &n, &m, &mod);
  printf("%lld\n", exlucas(n, m));
  return 0;
}