题目描述
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.
题意概述
求${n \choose k} \bmod m$。
数据范围:$1 \le n \le 10^{18}, \; 0 \le k \le n, \; 2 \le m \le 10^6$。
算法分析
令
$$
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.
$$
由于模数两两互质,所以该方程组在模$m$意义下有唯一解。
考虑如何求$r_i$。实际上,我们要求的就是${n \choose k} \bmod p_i^{a_i}$。我们知道
$$
{n \choose k}={n! \over k!(n-k)!}
$$
那么只要求出$n! \bmod p_i^{a_i}, \; k! \bmod p_i^{a_i}, \; (n-k)! \bmod p_i^{a_i}$的值,就可以用逆元求出${n \choose k} \bmod p_i^{a_i}$。
对于如何求$n! \bmod p_i^{a_i}$,令
$$
f(n)=n! \bmod p_i^{a_i}
$$
由$x \equiv x+p_i^{a_i} \pmod{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}
$$
但是$k!, \; (n-k)!$在模$p_i^{a_i}$意义下可能不存在逆元,因此需要将$n!, \; k!, \; (n-k)!$中的$p_i$因子提取出来,求出逆元后再乘回去。
这样就得到了所有$r_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; }