【BZOJ2445】最大团【推公式】【中国剩余定理】【扩展Lucas】

【题目链接】

公式为：

设 ans = ∑(n! / ((d!)^(n/d)*(n/d)!))

则答案为m ^ ans

证明：

考虑现在有d * k个点，d代表每个团的点数，那么k就是个数了，记方案数为Ak。

然后现在又来了d个点，记方案数为Ak+1。（即现在有n = d * (k + 1)个点）

我们选择一个点，让这个点与其他的d - 1个点组成团，方案数为C(d * (k + 1) - 1, d - 1)，那么就剩下了k * d个点，方案数为Ak，那么有

Ak+1 = C(d * (k + 1) - 1, d - 1) * Ak

把组合数展开成阶乘形式，再给分子分母同时乘上d * (k + 1)，然后递推展开，用n = d * (k + 1)代换，就可以得到上面式子。

注意题目中的p是个质数，那么我们求出ans mod (p - 1)即可，但是p - 1不是质数，直接求逆元的方法行不通。

但是发现p - 1 = 2 * 13 * 5281 * 7283。

我们用类似扩展Lucas的方法，最后CRT合并就行了。

/* Telekinetic Forest Guard */
#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

typedef long long LL;

const int maxn = 1000005;

int cnt = 4, pr[] = {0, 2, 13, 5281, 7283};

inline int mul(int a, int b, int p) {
	return (LL)a * b % p;
}

inline int qpow(int a, int n, int p) {
	int res = 1;
	for(; n; n >>= 1, a = mul(a, a, p)) if(n & 1) res = mul(res, a, p);
	return res;
}

inline void exgcd(int a, int b, int &x, int &y) {
	b ? (exgcd(b, a % b, y, x), y -= a / b * x) : (x = 1, y = 0);
}

inline int inv(int a, int p) {
	int x, y;
	exgcd(a, p, x, y);
	if(x < 0) x += p;
	return x;
}

inline int fact(int n, int p) {
	if(n == 0 || n == 1) return 1;
	int res = 1;
	for(int i = 2; i <= n && i <= p; i++) if(i % p) res = mul(res, i, p);
	res = qpow(res, n / p, p);
	for(int k = n % p, i = 2; i <= n && i <= k; i++) if(i % p) res = mul(res, i, p);
	return mul(res, fact(n / p, p), p);
}

int tot, di[maxn];

inline void getdiv(int n) {
	tot = 0;
	for(int i = 1; i * i <= n; i++) if(n % i == 0)
		di[++tot] = i;
}

inline int calc(int n, int pi, int p) {
	int factn = fact(n, pi), k = 0, res = 0;
	for(int x = n / pi; x; x /= pi) k += x;
	for(int i = 1; i <= tot; i++) {
		int factd = fact(di[i], pi), factnd = fact(n / di[i], pi);

		int tmpk = 0;
		for(int x = di[i] / pi; x; x /= pi) tmpk += x;
		tmpk *= n / di[i];
		for(int x = n / di[i] / pi; x; x /= pi) tmpk += x;
		(res += mul(mul(mul(factn, inv(qpow(factd, n / di[i], pi), pi), pi), inv(factnd, pi), pi), qpow(pi, k - tmpk, pi), pi)) %= pi;

		if(di[i] * di[i] != n) {
			swap(factd, factnd);
			tmpk = 0;
			for(int x = n / di[i] / pi; x; x /= pi) tmpk += x;
			tmpk *= di[i];
			for(int x = di[i] / pi; x; x /= pi) tmpk += x;
			(res += mul(mul(mul(factn, inv(qpow(factd, di[i], pi), pi), pi), inv(factnd, pi), pi), qpow(pi, k - tmpk, pi), pi)) %= pi;
		}
	}
	return mul(mul(res, p / pi, p), inv(p / pi, pi), p);
}

int main() {
	int T, p = 999999599;
	for(scanf("%d", &T); T; T--) {
		int n, m; scanf("%d%d", &n, &m);

		getdiv(n);

		int ans = 0;
		for(int i = 1; i <= 4; i++)
			(ans += calc(n, pr[i], p - 1)) %= p - 1;
		printf("%d\n", qpow(m, ans, p));
	}
	return 0;
}