将矩阵上每格对函数\(f\)值的贡献的系数画出来,会发现它是个斜着的正方形,中心系数等于\(s\),往外\(1\)层系数减\(1\),减到\(0\)为止,如下图。

CodeForces 263E Rhombus-LMLPHP

枚举中心点,可以得到\(\mathrm O(nm)\)个这样的正方形。现在问题是怎么快速求出所有正方形的系数乘格子内的数之和,即所有中心点的函数值。假设我们已经知道了\(f(i,j)\),现在想知道\(f(i,j+1)\)。不妨把这两个正方形画出来,看它们相差什么。

CodeForces 263E Rhombus-LMLPHP

如上图,是\(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\)为例,先画出相邻两个三角形。

CodeForces 263E Rhombus-LMLPHP

如上图,是两个相邻的\(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;
}
05-22 14:51