exLucas 学习笔记

发表于 更新于 分类于 阅读次数：

大质数取模有很多奇妙的姿势，其中比较菜的一个就是裸exLucas（然而我还是不会）。 exLucas适用于模数分解后不大的情况。 大概姿势是这样的：模数$P$分解为$p_1^{a_1} p_2^{a_2}…$，分别求解后用CRT合并。 那么求解$\binom{n+m}{n}$，只需要求解形如$n! \mod p^{a}$的式子。 令$n!=p^kc$，那么$k$很好求，只需要考虑$c$。注意到$1,2,3,…,p^a$在模意义下与$p^a+1,p^a+2,…,p^{2a}$完全一致，只需要预处理$p^a$之内，不是$p$倍数的积。递归即可。 例题：bzoj #3738. [Ontak2013]Kapitał 考虑分解成$2^a5^bc$，之后直接把上下的$10$消掉，套个exLucas上去就好了。代码参考（照抄）了Claris的实现，十分精妙。具体细节见代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include <bits/stdc++.h>

using std::cerr;
using std::endl;

#define int long long

int P = 1e9;

int exgcd(int a, int b, int &x, int &y) {
  if (!b) return x = 1, y = 0, a;
  int g = exgcd(b, a % b, y, x);
  return y -= a / b * x, g;
}

inline int get_inv(int val, int mod) {
  int x, y;
  exgcd(val, mod, x, y);
  x = (x % mod + mod) % mod;
  return x;
}

inline int fpow(int x, int y, int mod) {
  int ret = 1;
  for ( ; y; y >>= 1, x = x * x % mod)
    if (y & 1) ret = ret * x % mod;
  return ret;
}

int MOD;

// a * 2 ^ b
struct Node {
  int a, b;
  Node() { a = 1, b = 0; }
  Node(int x, int y) {
    a = x, b = y;
  }
  friend Node operator*(Node x, Node y) {
    return Node(x.a * y.a % MOD, x.b + y.b);
  }
  friend Node operator/(Node x, Node y) {
    return Node(x.a * get_inv(y.a, MOD) % MOD, x.b - y.b);
  }
} now[2];

int n, m, K, pw[1953125 + 1], del, res[2], ans, B;

Node calc(int val) {
  if (!val) return Node(1, 0);
  return Node(pw[val % MOD] * fpow(pw[MOD], val / MOD, MOD), val / B) * calc(val / B);
}

signed main() {
  std::cin >> n >> m >> K;
  MOD = 512, B = 2;
  for (int i = pw[0] = 1; i <= MOD; ++i)
    pw[i] = pw[i - 1] * (i % 2 ? i : 1) % MOD;
  pw[MOD] = pw[MOD - 1];
  now[0] = calc(n + m) / calc(n) / calc(m);
  MOD = 1953125, B = 5;
  for (int i = pw[0] = 1; i <= MOD; ++i)
    pw[i] = pw[i - 1] * (i % 5 ? i : 1) % MOD;
  pw[MOD] = pw[MOD - 1];
  now[1] = calc(n + m) / calc(n) / calc(m);
  del = std::min(now[0].b, now[1].b);
  while (del--) {
    MOD = 512, now[0] = now[0] / Node(5, 1);
    MOD = 1953125, now[1] = now[1] / Node(2, 1);
  }
  MOD = 512, res[0] = now[0].a * fpow(2, now[0].b, MOD) % MOD;
  MOD = 1953125, res[1] = now[1].a * fpow(5, now[1].b, MOD) % MOD;
  ans = (1953125 * get_inv(1953125, 512) % P * res[0] % P + 512 * get_inv(512, 1953125) % P * res[1] % P) % P;
  int T = 1;
  for (int i = 1; i <= K; ++i) T = T * 10;
  ans %= T;
  std::cout << std::setw(K) << std::setfill('0') << ans << endl;
  return 0;
}