题意

原题戳这里。

给你一个 (n imes m) 的矩阵 (A) ，一开始全是 (0) 。

然后你可以对这个矩阵进行 (k) 次操作：

每次选择一个宫格 ((i, j)) ，将第 (i) 行和第 (j) 列的状态翻转。（注意 ((i, j)) 这个点的状态翻转两次然后不会改变）

问最后使得 (A) 中恰好有 (S) 个宫格被染黑的方案数。

(1 le n, m, k le 3000)

题解

思路很简单，首先我们考虑最后如果有 (a) 行被染了奇数次， (b) 列被染了奇数次，那么最后被染黑（变成 (1) ）的格子数就是 (a * m + b * n - 2 * a * b) 。

我们可以暴力枚举 (a, b) 然后判断它得到面积是否为 (S) 。

然后我们就要考虑计算一个数，表示 (k) 次操作后，恰好有 (a) 行且 (b) 列被染黑的方案数。

我们很容易发现行和列的方案是互不影响的，然后我们可以考虑计算 (k) 次操作后，恰好有 (a) 行是染黑的方案数。

这个很显然可以用一个 (dp) 来计算，令 (f_{i, j}) 表示第 (i) 次操作后，恰好有 (j) 行染黑的方案数。

[f_{i,j} = f_{i - 1, j - 1} imes (n - j + 1) + f_{i - 1,j + 1} imes(j + 1) ]

这个转移意义十分的显然（就是忽略每行的区别，直接计算有几种方案使得它发生变化）

对与列的 (dp) 也是一样的。然后就能很轻松的做完啦qwq

复杂度是 (O(k(n + m) + nm)) 的。

总结

忽略一些表面上的区别，把本质相同的当做一个状态去转移就行了。

后记

我总认为这个题可以把 (n, m, k) 出到 (10^5) 乃至 (10^6) 。

对于 (a, b) 的枚举，我们可以单枚举一个，然后快速算出另外一个。（因为是一次函数，所以一一对应）

然后瓶颈就在求那个每行的方案数上了。

那个应该可以用生成函数等高端计数技巧来进行优化，目前我还没有找出规律qwq

其实可以用线性代数魔法优化到 (O(n log n log k)) ，就是特征多项式用 (NTT) 转移即可。

代码

#include <bits/stdc++.h>

#define For(i, l, r) for(register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for(register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << x << endl
#define DEBUG(...) fprintf(stderr, __VA_ARGS__)

using namespace std;

typedef long long ll;

inline bool chkmin(int &a, int b) {return b < a ? a = b, 1 : 0;}
inline bool chkmax(int &a, int b) {return b > a ? a = b, 1 : 0;}

inline int read() {
	int x = 0, fh = 1; char ch = getchar();
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') fh = -1;
	for (; isdigit(ch); ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48);
	return x * fh;
}

const int N = 3e3 + 1e2, Mod = 1e9 + 7;

void File() {
#ifdef zjp_shadow
	freopen ("binary-flips.in", "r", stdin);
	freopen ("binary-flips.out", "w", stdout);
#endif
}

ll dpn[N][N], dpm[N][N];

int main () {

	File();

	for (int cases = read(); cases; -- cases) {

		int n = read(), m = read(), k = read(), s = read();

		dpn[0][0] = 1;
		For (i, 1, k) { 
			For (j, 1, i)
				dpn[i][j] = (dpn[i - 1][j + 1] * (j + 1) + dpn[i - 1][j - 1] * (n - j + 1)) % Mod;
			dpn[i][0] = dpn[i - 1][1];
		}

		dpm[0][0] = 1;
		For (i, 1, k) { 
			For (j, 1, i) 
                dpm[i][j] = (dpm[i - 1][j + 1] * (j + 1) + dpm[i - 1][j - 1] * (m - j + 1)) % Mod;
			dpm[i][0] = dpm[i - 1][1];
		}

		int ans = 0;
		For (a, 0, k) if ((k & 1) == (a & 1)) For (b, 0, k) if ((k & 1) == (b & 1) && a * m + b * n - 2 * a * b == s) {
			ans = (ans + dpn[k][a] * dpm[k][b]) % Mod;
		}
		printf ("%d
", ans);

 }

	return 0;
}