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_{Tin 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)^2le O(n^2sum_{i=1}^n {1over 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("
")
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
", ans);
}