# Ceizenpok’s Formula

## 题目描述

Dr. Ceizenp’ok from planet i1c5l became famous across the whole Universe thanks to his recent discovery – the Ceizenpok’s formula. This formula has only three arguments: $n, k$ and $m$, and its value is a number of $k$-combinations of a set of $n$ modulo $m$.
While the whole Universe is trying to guess what the formula is useful for, we need to automate its calculation.

## 算法分析

$$x={n \choose k}, \; m=p_1^{a_1}p_2^{a_2} \cdots p_q^{a_q}$$

$$\left\{ \begin{array}{c} x \equiv r_1 \pmod{p_1^{a_1}} \\ x \equiv r_2 \pmod{p_2^{a_2}} \\ \cdots \\ x \equiv r_q \pmod{p_q^{a_q}} \end{array} \right.$$

$${n \choose k}={n! \over k!(n-k)!}$$

$$f(n)=n! \bmod p_i^{a_i}$$

$$f(n)=\left(f\left(\left\lfloor{n \over p_i}\right\rfloor\right)\cdot p_i^{\lfloor n/p_i \rfloor}\cdot\left(\prod_{i \in [1, p_i^{a_i}], \; p_i\not\mid i}i\right)^{\lfloor n/p_i^{a_i} \rfloor}\cdot\prod_{i \in [1, n \bmod p_i^{a_i}], \; p_i\not\mid i}i\right)\bmod p_i^{a_i}$$

## 代码

/*
* Your lucky number has been disconnected.
*/

#include <algorithm>
#include <cstdio>
#include <cstring>

#define int long long

int power(int a, int b, int m) {
int ret = 1;
for (; b; b >>= 1)
b & 1 && ((ret *= a) %= m), (a *= a) %= m;
return ret;
}

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

int get_inv(int a, int b) {
int x, y;
ex_gcd(a, b, x, y);
return (x + b) % b;
}

int get_fac(int n, int p, int k) {
int m = power(p, k, 1e9), ret = 1;
for (; n; n /= p) {
if (n / m) {
int rec = 1;
for (int i = 2; i < m; ++i)
if (i % p)
(rec *= i) %= m;
(ret *= power(rec, n / m, m)) %= m;
}
for (int i = n % m; i > 1; --i)
if (i % p)
(ret *= i) %= m;
}
return ret;
}

int calc(int N, int K, int M, int p, int k) {
int a = get_fac(N, p, k), b = get_fac(K, p, k), c = get_fac(N - K, p, k),
cnt = 0;
for (int i = N; i; i /= p)
cnt += i / p;
for (int i = K; i; i /= p)
cnt -= i / p;
for (int i = N - K; i; i /= p)
cnt -= i / p;
int m = power(p, k, 1e9),
ret = a * get_inv(b, m) % m * get_inv(c, m) % m * power(p, cnt, m) % m;
return ret * (M / m) % M * get_inv(M / m, m) % M;
}

signed main() {
int N, K, M, ans = 0;
scanf("%lld%lld%lld", &N, &K, &M);
for (int i = 2, t = M; t > 1; ++i)
if (!(t % i)) {
int k = 0;
for (; !(t % i); ++k, t /= i)
;
(ans += calc(N, K, M, i, k)) %= M;
}
printf("%lld\n", ans);
return 0;
}


418 I'm a teapot