【模板】扩展卢卡斯定理/exLucas

题目背景

这是一道模板题。

题目描述

求

\[{\mathrm{C}}_n^m \bmod{p} \]

其中 \(\mathrm{C}\) 为组合数。

输入格式

一行三个整数 \(n,m,p\) ，含义由题所述。

输出格式

一行一个整数，表示答案。

样例 #1

样例输入 #1

5 3 3

样例输出 #1

样例 #2

样例输入 #2

666 233 123456

样例输出 #2

提示

对于 \(100 \%\) 的数据，\(1 \le m \le n \le {10}^{18}\)，\(2 \le p \le {10}^6\)，不保证 \(p\) 是质数。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <cstdio>

typedef long long ll;

using namespace std;

ll n, m, sum, p;
ll a[5211314];

class Exlucas{
	long long c[5211314], r[5211314], d[5211314], tot;
	private:
		inline ll QuickPower(ll a, ll k, ll mod) {
			ll ans = 1, base = a;
			while (k != 0) {
				if (k & 1 == 1) {
					ans = ans * base % mod;
				}
				base = base * base % mod;
				k >>= 1;
			}
			return ans;
		}
		inline ll Exgcd(ll a, ll b, ll &x, ll &y) {
			if (b == 0) {
				x = 1, y = 0;
				return a;
			}
			ll x1, y1, gcd;
			gcd = Exgcd(b, a % b, x1, y1);
			x = y1;
			y = x1 - a / b * y1;
			return gcd;
		}
		inline ll Inv(ll a, ll mod) {
			ll x = 0, y = 0;
			Exgcd(a, mod, x, y);
			return (x % mod + mod) % mod;
		}
		inline ll Fac(ll n, ll p, ll pk) {
			if (n == 0) return 1;
			ll ans = 1;
			for (ll i = 1; i <= pk; ++ i) {
				if (i % p != 0) {
					ans = ans * i % pk;
				}
			}
			ans = QuickPower(ans, n / pk, pk);
			for (ll i = 1; i <= n % pk; ++ i) {
				if (i % p != 0) {
					ans = ans * i % pk;
				}
			}
			return ans * Fac(n / p, p, pk) % pk;
		}
		inline ll GetC(ll n, ll m, ll p, ll pk) {
			if (n < m) return 0;
			ll x = Fac(n, p, pk), y = Fac(m, p, pk), z = Fac(n - m, p, pk);
			ll k = n - m, cnt = 0;
			while (n != 0) n /= p, cnt += n;
			while (m != 0) m /= p, cnt -= m;
			while (k != 0) k /= p, cnt -= k;
			return x * Inv(y, pk) % pk * Inv(z, pk) % pk * QuickPower(p, cnt, pk) % pk;
		}
		inline ll CRT() {
			ll M = 1, ans = 0;
			for (ll i = 1; i <= tot; ++ i) {
				M *= c[i];
			}
			for (ll i = 1; i <= tot; ++ i) {
				d[i] = M / c[i];
			}
			for (ll i = 1; i <= tot; ++ i) {
				ans = (ans + r[i] * d[i] % M * Inv(d[i], c[i]) % M )% M;
			}
			return (ans % M + M) % M;
		}
	public:
		inline ll exlucas(ll n, ll m, ll p) {
			tot = 0;
			ll tmp = sqrt(p);
			for (ll i = 2; i <= tmp && p >= 1; ++ i) {
				ll pk = 1;
				while (p % i == 0) {
					p /= i, pk *= i;
				}
				if (pk > 1) {
					tot ++;
					c[tot] = pk;
					r[tot] = GetC(n, m, i, pk);
				}
			}
			if (p > 1) {
				tot ++;
				c[tot] = p;
				r[tot] = GetC(n, m, p, p);
			}
			cout << CRT() << endl;
			return CRT();
		}
}ask; 

int main() {
	cin >> n >> m >> p;
	ask.exlucas(n, m, p);
	return 0;
}