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的实现,十分精妙。具体细节见代码。

#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;
}

发表评论

电子邮件地址不会被公开。必填项已用 * 标注