@codechef

目录

• @desription@
• @solution@
• @accepted code@
• @details@

---

@desription@

给定一个 R * C 表示高度的矩阵 A，另一个 H * W 的矩阵 B，与一个坐标 (x, y)。

这个 R * C 的矩阵中的一个 H * W 的子矩阵，记这个子矩阵中某一个格子的差异值为 = (该方格相对于该子矩阵中的 (x, y) 的高度 - B 中对应方格的数值差)^2。
该子矩阵的差异值为所有格子的差异值之和。

求差异值前 K 小的子矩阵。输出其左上角坐标。

input
第一行输入两个整数 R，C（1 <= R, C <= 666）。
接下来 R 行每行 C 个数字，描述了矩阵 A。
接下来输入两个整数 H，W（1 <= H <= (R+1)/2，1 <= W <= (C+1)/2）
接下来 H 行每行 W 个数字，描述了矩阵 B。
接下来输入一个坐标 (x, y)。
最后输入一个整数 K（1 <= K <= 1000）。
保证所有矩阵元素的绝对值 <= 30。

output
输出 K 行，每行三个整数 (br, bc, s)，表示你选择的子矩阵左上角为 (br, bc)，其差异值为 s。
如果 s 相同，取 br 较小的；如果 br 相同，取 bc 较小的。

sample input
5 8
6 1 8 6 0 0 8 1
6 5 8 6 8 6 8 1
6 5 8 6 8 6 8 2
6 5 8 6 8 6 8 2
6 5 8 6 8 6 0 2
3 4
-1 2 2 2
-1 2 0 2
-1 2 2 2
2 3
4
sample output
2 2 8
3 2 8
2 4 11
3 4 75

@solution@

不妨去求解每一个左上角对应的差异值。我们考虑以 (i, j) 作为左上角时，差异值的表达式。

[sum_{p=1}^{H}sum_{q=1}^{W}(a[i+p-1][j+q-1] - a[i+x-1][j+y-1] - b[p][q])^2 ]

拆开可可以得到，它是由 6 部分组成：
第一部分：(H*W*a[i+x-1][j+y-1]^2)。
第二部分：(-2*a[i+x-1][j+y-1]*sum_{p=1}^{H}sum_{q=1}^{W}a[i+p-1][j+q-1])
第三部分：(2*a[i+x-1][j+y-1]*sum_{p=1}^{H}sum_{q=1}^{W}b[p][q])
第四部分：(sum_{p=1}^{H}sum_{q=1}^{W}a[i+p-1][j+q-1]^2)
第五部分：(sum_{p=1}^{H}sum_{q=1}^{W}b[p][q]^2)
第六部分：(-2*sum_{p=1}^{H}sum_{q=1}^{W}a[i+p-1][j+q-1]*b[p][q])

第 1 个可以 O(1) 算。第 2, 3, 4, 5 个可以做个前缀和然后 O(1) 算。
考虑第 6 个。

令 (c[i][j] = -2*sum_{p=1}^{H}sum_{q=1}^{W}a[i+p-1][j+q-1]*b[p][q])。
可以发现它非常的像卷积，但是不是 a 的下标与 b 的下标相加等于 c 的下标，而是 c 的下标与 b 的下标相加得到 a 的下标 - 1。

通过对 a 的翻转，移位，并在最后对 c 进行翻转，得到一个转换过后的问题：

[c[i+p][j+q] = sum_{p=1}^{H}sum_{q=1}^{W}a[i][j]*b[p][q] ]

令多项式 (A = sum a[i][j]*x^i*y^j)，同理得到 (B, C)。则 (C = A*B)。
这个问题虽然是卷积，但是有两个变元。不能用一般的 FFT 直接解决。

两种等价的方法，一个是插值时插若干二元组 (wi, wj)，将 wi 相同的放在一起对 wj 插值，过后再将 wj 相同的放在一起对 wi 插值。逆变换同理。
这种方法具有一定的拓展性，因为它将一个两变元多项式完全转换为了点值表示。

第二个要直观一些：
我们令 (A[i] = sum a[i][j]*x^j)，同理得到 (B[i], C[i])。
则有：(C[i+p] = sum A[i]*B[p])。

对 (A[i], B[i], C[i]) 作 FFT，得到 (C[i+p][j] = sum A[i][j]*B[p][j])。
交换横纵坐标的位置，得到 (C[j][i+p] = sum A[j][i]*B[j][p])。
这就是一个经典的卷积问题。

总时间复杂度 (O(n^2log n))，其中 n 与 R,C 同阶。

@accepted code@

#include<cstdio>
#include<algorithm>
using namespace std;
const int G = 3;
const int MAXN = 666;
const int MAXM = 2048 + 5;
const int MOD = 998244353;
int pow_mod(int b, int p) {
	int ret = 1;
	while( p ) {
		if( p & 1 ) ret = 1LL*ret*b%MOD;
		b = 1LL*b*b%MOD;
		p >>= 1;
	}
	return ret;
}
int pw[15], ipw[15];
void init() {
	for(int i=1;i<15;i++)
		pw[i] = pow_mod(G, (MOD - 1)/(1<<i)), ipw[i] = pow_mod(pw[i], MOD - 2);
}
void ntt(int *A, int n, int type) {
	for(int i=0,j=0;i<n;i++) {
		if( i < j ) swap(A[i], A[j]);
		for(int l=(n>>1);(j^=l)<l;l>>=1);
	}
	for(int k=1;(1<<k)<=n;k++) {
		int s = (1<<k), t = (s>>1);
		int u = (type == 1) ? pw[k] : ipw[k];
		for(int i=0;i<n;i+=s)
			for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {
				int x = A[i + j], y = 1LL*p*A[i + j + t]%MOD;
				A[i + j] = (x + y)%MOD, A[i + j + t] = (x + MOD - y)%MOD;
			}
	}
	if( type == -1 ) {
		int inv = pow_mod(n, MOD - 2);
		for(int i=0;i<n;i++)
			A[i] = 1LL*A[i]*inv%MOD;
	}
}
struct node{int x, y, ans; node(int _x=0, int _y=0, int _a=0):x(_x), y(_y), ans(_a){} }arr[MAXN*MAXN];
bool operator < (node a, node b) {
	if( a.ans == b.ans ) return a.x == b.x ? a.y < b.y : a.x < b.x;
	else return a.ans < b.ans;
}
int A[MAXM][MAXM], B[MAXM][MAXM], ans[MAXM][MAXM];
int M[MAXN][MAXN], s1[MAXN][MAXN], s2[MAXN][MAXN], S1, S2;
int getsum1(int x1, int y1, int x2, int y2) {
	int ret = s1[x2][y2];
	if( x1 > 0 ) ret -= s1[x1-1][y2];
	if( y1 > 0 ) ret -= s1[x2][y1-1];
	if( x1 > 0 && y1 > 0 ) ret += s1[x1-1][y1-1];
	return ret;
}
int getsum2(int x1, int y1, int x2, int y2) {
	int ret = s2[x2][y2];
	if( x1 > 0 ) ret -= s2[x1-1][y2];
	if( y1 > 0 ) ret -= s2[x2][y1-1];
	if( x1 > 0 && y1 > 0 ) ret += s2[x1-1][y1-1];
	return ret;
}
int R, C, H, W, K, x, y;
int main() {
	init();
	scanf("%d%d", &R, &C);
	for(int i=0;i<R;i++)
		for(int j=0;j<C;j++) {
			scanf("%d", &M[i][j]), s1[i][j] = M[i][j], s2[i][j] = M[i][j]*M[i][j];
			if( i > 0 ) s1[i][j] += s1[i-1][j], s2[i][j] += s2[i-1][j];
			if( j > 0 ) s1[i][j] += s1[i][j-1], s2[i][j] += s2[i][j-1];
			if( i > 0 && j > 0 ) s1[i][j] -= s1[i-1][j-1], s2[i][j] -= s2[i-1][j-1];
			A[R-i-1][C-j-1] = M[i][j];
		}
	scanf("%d%d", &H, &W);
	for(int i=0;i<H;i++)
		for(int j=0;j<W;j++)
			scanf("%d", &B[i][j]), S1 += B[i][j], S2 += B[i][j]*B[i][j], B[i][j] = -2*B[i][j];
	scanf("%d%d", &x, &y); x--, y--;
	int len; for(len = 1;len <= 2*R || len <= 2*C;len <<= 1);
	for(int i=0;i<R;i++) ntt(A[i], len, 1);
	for(int i=0;i<H;i++) ntt(B[i], len, 1);
	for(int i=0;i<len;i++)
		for(int j=0;j<i;j++)
			swap(A[i][j], A[j][i]), swap(B[i][j], B[j][i]);
	for(int i=0;i<len;i++) ntt(A[i], len, 1), ntt(B[i], len, 1);
	for(int i=0;i<len;i++)
		for(int j=0;j<len;j++)
			ans[j][i] = 1LL*A[j][i]*B[j][i]%MOD;
	for(int i=0;i<len;i++) ntt(ans[i], len, -1);
	for(int i=0;i<len;i++)
		for(int j=0;j<i;j++)
			swap(ans[i][j], ans[j][i]);
	for(int i=0;i<R;i++) ntt(ans[i], len, -1);
	for(int i=0;i<=R-H;i++)
		for(int j=0;j<=C-W;j++) {
			ans[R-i-1][C-j-1] += S2;
			ans[R-i-1][C-j-1] += getsum2(i, j, i+H-1, j+W-1);
			ans[R-i-1][C-j-1] += M[i+x][j+y]*M[i+x][j+y]*H*W;
			ans[R-i-1][C-j-1] += 2*M[i+x][j+y]*S1;
			ans[R-i-1][C-j-1] -= 2*M[i+x][j+y]*getsum1(i, j, i+H-1, j+W-1);
		}
	int cnt = 0;
	for(int i=0;i<=R-H;i++)
		for(int j=0;j<=C-W;j++)
			arr[cnt++] = node(i + 1, j + 1, ans[R-i-1][C-j-1]%MOD);
	sort(arr, arr + cnt);
	scanf("%d", &K);
	for(int i=0;i<K;i++)
		printf("%d %d %d
", arr[i].x, arr[i].y, arr[i].ans);
}

@details@

双变元什么的……以前还真的没遇到过。
不过如果按第二种理解方法的话，也不过是一些代数变形的 trick 罢了。

本题还有一种常数较大的方法：把 (i, j) 这一个坐标视为一个 C 进制数 i + j*C，其中 i 是第 1 位，j 是第 2 位。然后再做单变元的卷积。