拓展卢卡斯定理（伪）

题目

求 (a^{nchoose m}pmod{999999599})

数据范围均为1e9

思路

转化

首先观察到999999599是一个质数。根据费马小定理，我们只需要先求出 ({nchoose m}mod{999999598})，然后使用快速幂算法即可求得答案。

卢卡斯定理？

如果你还不知道卢卡斯定理是什么，可以阅读我的上一篇博文。

看到组合数取模这样的式子，自然想到用卢卡斯定理来求解。但是卢卡斯定理有一个限制：模数必须是质数。

拓展卢卡斯定理？

拓展卢卡斯定理是一个对卢卡斯定理的推广，可以解决模数不是质数的情况。关于它的详细资料可以参照卢卡斯定理 - OI Wiki。

正解

然而这题不用拓展卢卡斯定理就能做啊Σ(っ °Д °;)っ或者说这是一个拓展卢卡斯定理的简化情况。

我们对999999598进行质因数分解：

嗯？正好四个质因数？那么我们可以使用如下方法求解：

[egin{cases}Cequiv{nchoose m}pmod{2}\Cequiv{nchoose m}pmod{13}\Cequiv{nchoose m}pmod{5281}\Cequiv{nchoose m}pmod{7283}end{cases} ]

也就是说，我们可以先求出({nchoose m}mod{2})，({nchoose m}mod{13})，({nchoose m}mod{5281})和({nchoose m}mod{7283})，然后使用中国剩余定理来解决！

代码

注：代码中那些莫名其妙的常数是中国剩余定理的中间计算结果

#define proname "silk"
#include <cstdio>

#define ll long long
#define re register
#define il inline
#define gc getchar
#define pc putchar

template <class T>
void read(T &x) {
  re bool f = 0;
  re char c = gc();
  while ((c < '0' || c > '9') && c != '-') c = gc();
  if (c == '-') f = 1, c = gc();
  x = 0;
  while (c >= '0' && c <= '9') x = x * 10 + (c ^ 48), c = gc();
  f && (x = -x);
}

template <class T>
void print(T x) {
  if (x < 0) pc('-'), x = -x;
  if (x >= 10) print(x / 10);
  pc((x % 10) ^ 48);
}

template <class T>
void prisp(T x) {
  print(x);
  pc(' ');
}
template <class T>
void priln(T x) {
  print(x);
  pc('
');
}

ll a, n, m;

ll fac[10000];

ll pow(ll b, ll t, int p) {
  ll r = 1;
  for(; t; t>>=1, b = (b * b) % p) if (t & 1) r = (r * b) % p;
  return r;
}

ll C(ll n, ll k, int p) {
  if (k > n) return 0;
  return (fac[n] * pow(fac[k], p - 2, p) % p) * pow(fac[n - k], p - 2, p) % p;
}

ll lucas(ll n, ll k, ll p) {
  if (!k) return 1;
  return C(n % p, k % p, p) * lucas(n / p, k / p, p) % p;
}

ll calc(ll n, ll k, int p) {
  fac[0] = 1;
  for (int i = 1; i < p; ++i) fac[i] = (fac[i - 1] * i) % p;
  return lucas(n, k, p);
}

const int P = 999999599, M = P - 1;

int main() {
  // freopen(proname".in", "r", stdin);
  // freopen(proname".out", "w", stdout);
  read(a);
  read(n);
  read(m);
  ll c1 = calc(n, m, 2);
  ll c2 = calc(n, m, 13);
  ll c3 = calc(n, m, 5281);
  ll c4 = calc(n, m, 7283);
  ll c = M / 2 * c1 % M;
  c += M / 13 * 8 * c2 % M;
  c %= M;
  c += M / 5281 * 4647 * c3 % M;
  c %= M;
  c += M / 7283 * 34 * c4 % M;
  c %= M;
  ll ans = pow(a, c, P);
  priln(ans);
}