USACO 2020 Open Contest, Platinum（exercise）（min-max容斥+计数dp）

https://gmoj.net/senior/#main/show/6545

考场上想到了min-max容斥，结果后面的dp用了复杂度劣的，多了个log，就只有75p了。

求所有lcm的乘积，转换为枚举一个质数\(p\)，求p出现的指数\(mod ~ mo-1\)

那么\(=\sum_{一种划分S} S含有的p的指数的max*f(S)\)
S这种划分就是把n分成若干环，环长和\(=n\)，\(f(S)\)表示分成这些环的方案数。

暴力：
做动态规划\(f[i][j]\)表示现在的指数max是\(i\)，环长和是\(j\)的方案数。

转移：从小到大枚举环长，再枚举这个环长的个数，乘上些系数，复杂度\(O(log_p(n)*n^2*log~n)\)
转移2：枚举最小元素所在的环长，乘上些系数，复杂度\(O(log_p(n)*n^2)\)

当然，这个dp的第一维也可以不要，枚举个\(k\)，求指数max至多是\(k\)的方案数，发现可以按环长含有p的指数从小到大添加新的环长，这样可以去掉\(log_p(n)\)的复杂度。

求\(max\)一类问题往往不好做，利用min-max容斥，
变为\(\sum_{T}min(T)*(-1)^{|T|+1}*\sum_{T\in S}f(S)\)
\(=\sum_{T}min(T)*(-1)^{|T|+1}*f(T)*((S/T)的元素和)!\)

这时,\(min(T)\)显然要\(>0\)才有意义。
还是枚举个k，表示求\(min>=k\)的方案数，那么把\(p^k\)的倍数作为环长丢进dp，注意因为都是\(p^k\)的倍数，可以假设它们都除以了\(p^k\)，背包上界变为\(n/p^k\)。

总复杂度：\(O(\sum_{p是质数}~~~\sum_{k>0~and~p^k<=n}(n/p^k)^2\le O(n^2\sum_{i=1}^n {1\over i^2})\)

预处理组合数就可以处理模数不是质数的情况。

Code

#pragma GCC optimize(2)
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;

int gcd(int x, int y) {
	return !y ? x : gcd(y, x % y);
}

const int N = 7505;

int n, mo, mo2;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

int bp[N], p[N], p0;

ll c[N][N], fac[N];

ll f[N];

void dp2(int p) {
	fo(i, 0, n / p) f[i] = 0;
	f[0] = -1;
	fo(i, 1, n / p) {
		fo(j, 1, i) {
			f[i] = (f[i] - f[i - j] * c[i * p - 1][j * p - 1] % mo2 * fac[j * p - 1]) % mo2;
		}
	}
}

int main() {
	freopen("exercise.in", "r", stdin);
	freopen("exercise.out", "w", stdout);
	scanf("%d %d", &n, &mo);
	mo2 = mo - 1;
	fo(i, 0, n) {
		c[i][0] = 1;
		fo(j, 1, i) c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mo2;
	}
	fo(i, 2, n) fo(j, 2, n / i) bp[i * j] = 1;
	fo(i, 2, n) if(!bp[i]) p[++ p0] = i;
	fac[0] = 1; fo(i ,1, n) fac[i] = fac[i - 1] * i % mo2;
	ll ans = 1;
	fo(i, 1, p0) {
		int mx = log2(n) / log2(p[i]);
		ll sum = 0;
		ll sp = 1;
		fo(j, 1, mx) {
			sp *= p[i];
			dp2(sp);
			fo(k, 1, n / sp) sum = (sum + f[k] * fac[n - k * sp] % mo2 * c[n][k * sp]) % mo2;
		}
		if(sum < 0) sum += mo2;
		ans = ans * ksm(p[i], sum) % mo;
	}
	pp("%lld\n", ans);
}