将矩阵上每格对函数\(f\)值的贡献的系数画出来，会发现它是个斜着的正方形，中心系数等于\(s\)，往外\(1\)层系数减\(1\)，减到\(0\)为止，如下图。
枚举中心点，可以得到\(\mathrm O(nm)\)个这样的正方形。现在问题是怎么快速求出所有正方形的系数乘格子内的数之和，即所有中心点的函数值。假设我们已经知道了\(f(i,j)\)，现在想知道\(f(i,j+1)\)。不妨把这两个正方形画出来，看它们相差什么。
如上图，是\(2\)个相邻的\(s=3\)的正方形，紫色的数字是此格子内红色系数减蓝色系数。不难发现，左边一半是一个差都是\(-1\)、高为\(s=3\)、直角顶点在底边左侧的等腰直角三角形，右边一半是一个差都是\(1\)……底边右侧……。由此可以归纳出，设\(trl_{i,j}\)表示底边中点为\((i,j)\)、高为\(s\)、直角顶点在底边左侧的等腰直角三角形内元素之和，\(trr_{i,j}\)表示……底边右侧……，则\(f(i,j)=f(i,j-1)-trl_{i,j-1}+trr_{i,j}\)。假设我们已经知道了\(trl,trr\)数组，那么可以\(\forall i\in[s,n-s+1]\)，暴力用对角线前缀和\(\mathrm O(s)\)求出\(f(i,s)\)，然后\(\forall i\in[s,n-s+1],\forall j\in(s,m-s+1]\)，用上面那个关系式递推求出\(f(i,j)\)，总复杂度为\(\mathrm O(nm)\)。
现在问题转化为怎么快速求出\(trl,trr\)数组。以\(trl\)为例，先画出相邻两个三角形。
如上图，是两个相邻的\(s=3\)的三角形，紫色的数字是此格子内红色系数减蓝色系数。又不难发现，左边的三角形轮廓的差都是\(-1\)，右边一列的差都是\(1\)。这个轮廓的和和列的和可以维护列前缀和\(Sum1\)、副对角线前缀和\(Sum2\)（在同一条副对角线上当且仅当行列和相等）、主对角线前缀和\(Sum3\)（在同一条主对角线上当且仅当行列差相等）搞定。于是就有了一个\(trl_{i,j},trl_{i,j-1}\)的关系式。那么可以\(\forall i\in[s,n-s+1]\)，暴力用列前缀和\(\mathrm O(s)\)求出\(trl_{i,s}\)，然后\(\forall i\in[s,n-s+1],\forall j\in(s,m]\)，用关系式递推求出\(trl_{i,j}\)，总复杂度为\(\mathrm O(nm)\)。\(trr\)求法类似。
\(\mathrm O(nm)\)与处理完\(3\)个前缀和，\(trl,trr\)就可以\(\mathrm O(nm)\)求出来了，那么\(f\)也可以\(\mathrm O(nm)\)求出来了，最后找最大值即可。总复杂度\(\mathrm O(nm)\)。
下面贴代码：
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define mp make_pair
#define X first
#define Y second
const int N=1000,M=1000;
int n/*矩阵行数*/,m/*矩阵列数*/,s/*常数*/;
int a[N+1][M+1];//矩阵
int Sum1[M+1][N+1]/*列前缀和*/,Sum2[N+M+1][M+1]/*副对角线前缀和*/,Sum3[2*N+1][M+1]/*主对角线前缀和，由于行列差可能是负数，所以平移max(n,m)个单位*/;
int sum1(int x,int l,int r){return l>r?0:Sum1[x][r]-Sum1[x][l-1];}//列区间和
int sum2(int x,int l,int r){return l>r?0:Sum2[x][r]-Sum2[x][l-1];}//副对角线区间和
int sum3(int x,int l,int r){return l>r?0:Sum3[x+max(n,m)][r]-Sum3[x+max(n,m)][l-1];}//主对角线区间和
int trl[N+1][M+1],trr[N+1][M+1];//朝左、朝右三角形
int rhm[N+1][M+1];//正方形，即f函数
signed main(){
cin>>n>>m>>s;
for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)scanf("%lld",a[i]+j);
//预处理前缀和开始
for(int i=1;i<=m;i++)for(int j=1;j<=n;j++)Sum1[i][j]=Sum1[i][j-1]+a[j][i];
for(int i=1;i<=n+m;i++)for(int j=1;j<=n;j++)Sum2[i][j]=Sum2[i][j-1]+(1<=i-j&&i-j<=m?a[j][i-j]:0);
for(int i=-max(n,m);i<=max(n,m);i++)for(int j=1;j<=n;j++)Sum3[i+max(n,m)][j]=Sum3[i+max(n,m)][j-1]+(1<=j+i&&j+i<=m?a[j][j+i]:0);
//预处理前缀和结束
for(int i=s;i<=n-s+1;i++){//算trl
for(int j=1;j<=s;j++)trl[i][s]+=sum1(j,i-j+1,i+j-1);//暴力算边上的trl
for(int j=s+1;j<=m;j++)trl[i][j]=trl[i][j-1]+sum1(j,i-s+1,i+s-1)-sum2(i+j-s,i-s+1,i)-sum3(j-s-i,i+1,i+s-1);//用关系式递推其他trl
}
for(int i=s;i<=n-s+1;i++){//算trr，与trl类似
for(int j=m;j>=m-s+1;j--)trr[i][m-s+1]+=sum1(j,i-(m-j+1)+1,i+(m-j+1)-1);
for(int j=m-s;j;j--)trr[i][j]=trr[i][j+1]+sum1(j,i-s+1,i+s-1)-sum3(j+s-i,i-s+1,i)-sum2(i+j+s,i+1,i+s-1);
}
// for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)printf("tr(%lld,%lld)=(%lld,%lld)\n",i,j,trl[i][j],trr[i][j]);
for(int i=s;i<=n-s+1;i++){//算rhm
for(int j=0;j<s;j++)rhm[i][s]+=(s-j)*(sum2(i-j+s,i-j,i)+sum3(s-(i+j),i+1,i+j)+sum2(i+j+s,i,i+j-1)+sum3(s-(i-j),i-j+1,i-1));//暴力算边上的rhm
for(int j=s+1;j<=m-s+1;j++)rhm[i][j]=rhm[i][j-1]+trr[i][j]-trl[i][j-1];//用关系式递推其他rhm
}
// for(int i=s;i<=n-s+1;i++)for(int j=s;j<=m-s+1;j++)printf("rhm[%lld][%lld]=%lld\n",i,j,rhm[i][j]);
pair<int,pair<int,int> > mx(0,mp(0,0));
for(int i=s;i<=n-s+1;i++)for(int j=s;j<=m-s+1;j++)mx=max(mx,mp(rhm[i][j],mp(i,j)));//找最大值
cout<<mx.Y.X<<" "<<mx.Y.Y;
return 0;
}