@codechef - BUYLAND@ Buying Land
@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) 作为左上角时，差异值的表达式。
拆开可可以得到，它是由 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 进行翻转，得到一个转换过后的问题：
令多项式 \(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^2\log 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\n", arr[i].x, arr[i].y, arr[i].ans);
}
@details@
双变元什么的……以前还真的没遇到过。
不过如果按第二种理解方法的话，也不过是一些代数变形的 trick 罢了。
本题还有一种常数较大的方法：把 (i, j) 这一个坐标视为一个 C 进制数 i + j*C，其中 i 是第 1 位，j 是第 2 位。然后再做单变元的卷积。
