BZOJ4673 置换（动态规划）

BZOJ4673 置换（动态规划）

题目大意

对于由n个元素的置换P，我们定义(f(P))为最小的正整数k，使得(P^k)为恒等置换.

举个例子，当n=3时，考虑置换(P=(3 1 2))，那么(P^2=(2 3 1)，P^3=(1 2 3))，因此(f(P)=3).

你需要求出对于所有的n元素置换P，(f^2(P))的平均值是多少？

数据范围

$ n le 200, 10000 < p < 10^9 $

解题思路

老套路题了

一个显然的结论是$ f(P)=lcm(循环节大小)$

第一步

考虑 dp，(dp[x][y]) 表示大小为 x 的一个置换，(P(x)=y) 的方案数是多少

转移中枚举加入 a 个大小为 s 的循环节，组合数转移，分配标号要用可重排列，大小为 s 的循环节有 ((s-1)!) 个，即圆排列，同时要除以 (a!) 消除 s 间顺序的影响

[dp[x+a*s][lcm(y, s)] = dp[x][y]*prod_{i=1}^a{x + i*s choose s}*((s-1)!)^a * (a!)^{-1} ]

一点小优化

可重排列最后的形式为 (frac{n!}{prod x!})，因为我们每个 x 又乘上了 ((x-1)!)，最后答案除去了 (n!)，可化为 (frac{1}{prod x})，所以转移时直接乘上 (x^a) 的逆元和 (a!) 的逆元即可

正解

发现 lcm 可能很大，数组开都开不下，爆搜将所有可能的 lcm 都存下来，发现大概有 80 多万，仍然很 shi，可以用上我们的套路，发现大于 (sqrt n) 的质因数次数只能是一，把小于 (sqrt n) 的质数拼成的 lcm 存下来发现只有 2000 左右可以做了

我们将添加大小为 x 的置换的顺序按照 x 的大于 (sqrt n) 的质数排序，同层转移，类似 NOI2015寿司晚宴。同层之间肯定使 lcm 乘上这个质数，设 s 为大小为 x，lcm 为 y（已压缩且只含小于 (sqrt n) 的质数）且必须使用此层的置换的方案数，将 s 乘上 y * y 存到 (dp[x][y]) 中去，也就是将其看作 s * y * y 种方案数

带码

const int N = 2500;
int mx = 200;
ll res[N], tot, cnt;
int prime[] = {0, 2, 3, 5, 7, 11, 13, 17, 19};

void dfs(int x, int sum, ll mul) {
	if (sum > mx) return;
	if (x > 6) return res[++tot] = mul, void();
	dfs(x + 1, sum, mul);
	ll tmp = 1; 
	for (int i = 1;i <= mx; i++) {
		tmp *= prime[x];
		if (sum + tmp > mx) break;
		dfs(x + 1, sum + tmp, mul * tmp); 
	}
}

ll f[505][N], g[505][N], n, p;

inline ll lcm(ll a, ll b) {
	return a * b / __gcd(a, b);
}

inline int get(ll k, ll t) {
	t = lower_bound(res + 1, res + tot + 1, lcm(k, t)) - res;
	return t;
}

ll inv[N], rk[N], pp[N], P;
int main() {
	
	read(n), read(p); mx = n, P = p;
	dfs(1, 0, 1), sort(res + 1, res + tot + 1);
	inv[0] = inv[1] = 1; 
	for (int i = 2;i <= n; i++) inv[i] = (P - P / i) * inv[P % i] % P;
	for (int i = 1;i <= n; i++) {
		rk[i] = pp[i] = i;
		for (int j = 1;j <= 6; j++)
			while (pp[i] % prime[j] == 0) pp[i] /= prime[j];
	}
	sort(rk + 1, rk + n + 1, [](int a, int b) { return pp[a] < pp[b]; });
	
	
	
	ll ans = 0; f[0][1] = 1;
	for (int i = 1, j = 1; i <= n ; i = j) {
		for (int k = 0;k <= n; k++) 
			for (int q = 0;q <= tot; q++) 
				g[k][q] = f[k][q];
		while (pp[rk[i]] == pp[rk[j]]) {
			int nw = rk[j]; j++;
			for (int k = tot;k >= 1; k--) {
				int nxt = get(res[k], nw / pp[nw]);
				if (nxt > tot) continue;
//				printf ("%d %d %d
", res[k], nw / pp[nw], nxt);
				for (int a = n - nw; a >= 0; a--) {
					if (!g[a][k]) continue;
					int uc = 1, y = a + nw;
					ll div = inv[nw];
					for ( ; y <= n ;uc++, div = div * inv[nw] % p * inv[uc] % p, y += nw) 
						g[y][nxt] = (g[y][nxt] + g[a][k] * div) % p;
				}
			}
			
//			for (int i = 1;i <= tot; i++) 
//				printf ("%lld
", g[n][i]);
		}
		
		ll t = pp[rk[i]] * pp[rk[i]] % P;
		for (int k = 0;k <= n; k++)
			for (int q = 1;q <= tot; q++)
				f[k][q] = (f[k][q] + t * (g[k][q] - f[k][q])) % P;
	}
	
	for (int i = 1;i <= tot; i++) {
//		printf ("%lld
", f[n][i]);
		(ans += f[n][i] * res[i] % p * res[i]) %= p;
	}
	cout << (ans % P + P) % P << endl;
	
	return 0;
}

/*



*/