AtCoder Grand Contest 039 F: Min Product Sum

Posted 2020-11-28 pinkrabbit

题目传送门：AGC039F。

题意简述

有一个 (N imes M) 的矩阵，其中元素是 (1 sim K) 之间的整数。显然总共有 (K^{N M}) 种不同的矩阵。

定义一个矩阵的权值为：枚举每个元素 ((i, j))，将这个元素的所在行和所在列中的所有元素（共 (N + M - 1) 个）中的最小值求出来，记作 (f(i, j))，则权值为所有 (f(i, j)) 的乘积，即 (displaystyle prod_{substack{1 le i le N \ 1 le j le M}} f(i, j))。

求所有 (K^{N M}) 种不同的矩阵的权值之和，对 (mathrm{MOD}) 取模，其中 ({10}^8 le mathrm{MOD} le {10}^9)。

• (1 le N, M, K le 100)。

题解

我们将这个矩阵称为 (A)。

考虑枚举 (A) 的每一行的最小值和每一列的最小值，分别记作 (x_{1 sim N}) 和 (y_{1 sim M})。
也就是 (displaystyle x_i = min_{j = 1}^{M} A_{i, j}) 以及 (displaystyle y_j = min_{i = 1}^{N} A_{i, j})。

假设 (x, y) 已知，则迅速得到 (A) 的权值就为 (displaystyle prod_{substack{1 le i le N \ 1 le j le M}} min(x_i, y_j))。

我们考虑一种可以在已知 (x, y) 的条件下，在 (mathcal O (N + M)) 个乘法运算内求出此权值的方法：
将 (x, y) 从小到大排序，每次加入一个最小的元素，如果来自 (x) 则称为加入一行，否则称为加入一列。
假设已经加入 (i) 行 (j) 列，则加入一个元素 (t) 时产生的贡献为 (t^{M - j})（来自 (x)）或 (t^{N - i})（来自 (y)）。
将每次加入元素时产生的贡献全部相乘即可得到权值。

其原因是：因为从小到大加入，所以在第一次加入时本行/列就全部被这个元素限制了，直接乘上 (t) 的还未被限制的元素数量次方。

以上是从小到大考虑，还有一种从大到小考虑的方法：
每次加入一个最大的元素，假设已经加入 (i) 行 (j) 列，则加入一个元素 (t) 时产生的贡献为 (t^j)（来自 (x)）或 (t^i)（来自 (y)）。

原因相似，因为从大到小加入，所以每个元素是被行列较晚加入的那个限制的，这就导致影响到的元素数量为已加入的行/列个数。

---

接下来会用到这两种求法。以上是 (x, y) 固定的情况，现在我们考虑枚举 (x, y)。

首先可以发现 (x, y) 是有序的（不同顺序算不同种），但上述过程的 (x, y) 无序（不同顺序算同种，换句话说被人为排序了）。

这提示我们把 (x, y) 看作无序，从小到大加入元素，但是必须一次性加入所有一样的元素，以便用多重组合数转为有序的情况。

枚举 (x, y) 之后，其权值可以确定，但是还需要计算满足 (A) 的行列最小值恰好等于 (x, y) 的矩阵 (A) 的数量。

「恰好」提示我们需要容斥：假设不是恰好，而是 (A) 的行列最小值大于等于 (x, y)。

也就是已知 (x, y)，对于 (A_{i, j})，它的取值为 (max(x_i, y_j) sim K) 之间的正整数。

如果我们把权值在 ([1, K]) 内取反，可以发现「从小到大」枚举 (x, y) 中的元素就变成了「从大到小」，恰好对应上述第二种方法。

枚举 (c) 行 (d) 列容斥，则这些行列的限制就不是大于等于而是大于，或者说是大于等于「原限制加 (1)」，且容斥系数为 ({(-1)}^{c + d})。

可以把这 (c) 行 (d) 列混进从小到大加入 (x, y) 中的元素的过程中，也就是从小到大加入，并且分为是否「被容斥」两种类型。

考虑这样的 DP 过程：令 (mathrm{f}[t][i][j]) 表示加入了 (x, y) 中所有小于等于 (t) 的元素，所有可能情况下的系数之和。

注意必须一次加入所有相同元素以使用多重组合数确定还原成有序的 (x, y) 的方法数量。

先枚举 (a)，加入 (a) 行不被容斥的值为 (t) 的行：

[{color{red}{mathrm{f}[t][i + a][j]}} gets left[ {color{orange}{inom{N - i}{a}}} cdot {color{blue}{t^{a (M - j)}}} cdot {color{green}{{(K - t + 1)}^{a j}}} cdot {color{red}{mathrm{g}[i][j]}} ight] ]

然后枚举 (b)，加入 (b) 列不被容斥的值为 (t) 的列：

[{color{red}{mathrm{f}[t][i][j + b]}} gets left[ {color{orange}{inom{M - j}{b}}} cdot {color{blue}{t^{(M - i) b}}} cdot {color{green}{{(K - t + 1)}^{i b}}} cdot {color{red}{mathrm{g}[i][j]}} ight] ]

然后枚举 (c)，加入 (c) 行被容斥的值为 (t) 的行：

[{color{red}{mathrm{f}[t][i + c][j]}} gets left[ {color{orange}{inom{N - i}{c}}} cdot {color{blue}{t^{c (M - j)}}} cdot {color{green}{{(K - t)}^{c j}}} cdot {color{brown}{{(-1)}^c}} cdot {color{red}{mathrm{g}[i][j]}} ight] ]

最后枚举 (d)，加入 (d) 列被容斥的值为 (t) 的列：

[{color{red}{mathrm{f}[t][i][j + d]}} gets left[ {color{orange}{inom{M - j}{d}}} cdot {color{blue}{t^{(M - i) d}}} cdot {color{green}{{(K - t)}^{i d}}} cdot {color{brown}{{(-1)}^d}} cdot {color{red}{mathrm{g}[i][j]}} ight] ]

其中 (mathrm{g}) 就表示上一阶段的 DP 数组。可以滚动数组优化，具体实现详见代码，本题略有卡常。

下面是代码，时间复杂度为 (mathcal O (K N M (N + M)))：

#include <cstdio>
#include <algorithm>

typedef long long LL;
int Mod;
const int MN = 105, MK = 105;

int Binom[MN][MN], Pow[MK][MN];
int N, M, K;
int f[2][MN][MN];

int main() {
	scanf("%d%d%d%d", &N, &M, &K, &Mod);
	for (int i = 0; i <= std::max(N, M); ++i) {
		Binom[i][0] = 1;
		for (int j = 1; j <= i; ++j)
			Binom[i][j] = (Binom[i - 1][j - 1] + Binom[i - 1][j]) % Mod;
	}
	for (int i = 0; i <= K; ++i) {
		Pow[i][0] = 1;
		for (int j = 1; j <= std::max(N, M); ++j)
			Pow[i][j] = (LL)Pow[i][j - 1] * i % Mod;
	}
	int o = 0;
	f[0][0][0] = 1;
#define Z0(x) for (int i = 0; i <= N; ++i) for (int j = 0; j <= M; ++j) f[x][i][j] = 0;
	for (int t = 1; t <= K; ++t) {
		Z0(o ^ 1);
		for (int i = 0; i <= N; ++i) for (int j = 0; j <= M; ++j) if (f[o][i][j]) {
			int x = f[o][i][j], y = (LL)Pow[K - t + 1][j] * Pow[t][M - j] % Mod;
			for (int a = 0; i + a <= N; ++a, x = (LL)x * y % Mod)
				f[o ^ 1][i + a][j] = (f[o ^ 1][i + a][j] + (LL)Binom[N - i][a] * x) % Mod;
		} o ^= 1;
		Z0(o ^ 1);
		for (int i = 0; i <= N; ++i) for (int j = 0; j <= M; ++j) if (f[o][i][j]) {
			int x = f[o][i][j], y = (LL)Pow[K - t + 1][i] * Pow[t][N - i] % Mod;
			for (int b = 0; j + b <= M; ++b, x = (LL)x * y % Mod)
				f[o ^ 1][i][j + b] = (f[o ^ 1][i][j + b] + (LL)Binom[M - j][b] * x) % Mod;
		} o ^= 1;
		Z0(o ^ 1);
		for (int i = 0; i <= N; ++i) for (int j = 0; j <= M; ++j) if (f[o][i][j]) {
			int x = f[o][i][j], y = (LL)(Mod - Pow[K - t][j]) * Pow[t][M - j] % Mod;
			for (int c = 0; i + c <= N; ++c, x = (LL)x * y % Mod)
				f[o ^ 1][i + c][j] = (f[o ^ 1][i + c][j] + (LL)Binom[N - i][c] * x) % Mod;
		} o ^= 1;
		Z0(o ^ 1);
		for (int i = 0; i <= N; ++i) for (int j = 0; j <= M; ++j) if (f[o][i][j]) {
			int x = f[o][i][j], y = (LL)(Mod - Pow[K - t][i]) * Pow[t][N - i] % Mod;
			for (int d = 0; j + d <= M; ++d, x = (LL)x * y % Mod)
				f[o ^ 1][i][j + d] = (f[o ^ 1][i][j + d] + (LL)Binom[M - j][d] * x) % Mod;
		} o ^= 1;
	}
	printf("%d
", f[o][N][M]);
	return 0;
}

注意到编译器是不会使用 Barrett Reduction 算法优化对变量取模的时间常数的。
但是由于 AtCoder 实际上是支持 __uint128_t 的，所以可以使用手写的 Barrett Reduction 算法以大幅度减小常数。
官方题解中提到了一种每个 (t) 仅需分 (2) 个阶段进行转移的做法，也就是不需要使用容斥原理，应该也有一半的常数。