BZOJ4624 : 农场种植

设$A[i][j]=[a[i][j]=G]，B[i][j]=[b[i][j]=L]$，枚举右下角，则对应$(A-B)^2$的和就是匹配成功的格子数。

$(a-b)^2=a^2+b^2-2ab$，将矩阵展开成一维，然后将$B$翻转，用FFT求出$sum a imes b$即可。

时间复杂度$O(BRClog(RC))$。

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int N=530000,M=505;
int Case,T,r,c,n,m,i,j,k,pos[N],tmp,val,s[M][M],ans,sx,sy,c0,c1;char a[M][M],b[M][M];
struct comp{
  double r,i;comp(double _r=0,double _i=0){r=_r;i=_i;}
  comp operator+(const comp&x){return comp(r+x.r,i+x.i);}
  comp operator-(const comp&x){return comp(r-x.r,i-x.i);}
  comp operator*(const comp&x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}
  comp conj(){return comp(r,-i);}
}A[N],B[N];
const double pi=acos(-1.0);
void FFT(comp a[],int n,int t){
  for(int i=1;i<n;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
  for(int d=0;(1<<d)<n;d++){
    int m=1<<d,m2=m<<1;
    double o=pi*2/m2*t;comp _w(cos(o),sin(o));
    for(int i=0;i<n;i+=m2){
      comp w(1,0);
      for(int j=0;j<m;j++){
        comp&A=a[i+j+m],&B=a[i+j],t=w*A;
        A=B-t;B=B+t;w=w*_w;
      }
    }
  }
  if(t==-1)for(int i=0;i<n;i++)a[i].r/=n;
}
int main(){
  scanf("%d%d",&r,&c);
  for(i=0;i<r;i++)scanf("%s",a[i]);
  for(i=0;i<r;i++)for(j=0;j<c;j++){
    tmp=a[i][j]=='G';
    A[i*c+j]=comp(tmp,0);
    s[i][j]=tmp*tmp;
    if(i)s[i][j]+=s[i-1][j];
    if(j)s[i][j]+=s[i][j-1];
    if(i&&j)s[i][j]-=s[i-1][j-1];
  }
  for(k=1;k<r*c;k<<=1);k<<=1;
  j=__builtin_ctz(k)-1;
  for(i=0;i<k;i++)pos[i]=pos[i>>1]>>1|((i&1)<<j);
  FFT(A,k,1);
  scanf("%d",&Case);
  for(T=1;T<=Case;T++){
    scanf("%d%d",&n,&m);
    for(i=0;i<n;i++)scanf("%s",b[i]);
    for(i=0;i<k;i++)B[i]=comp();
    for(val=i=0;i<n;i++)for(j=0;j<m;j++){
      tmp=b[i][j]=='L';
      val+=tmp*tmp;
      B[(n-i-1)*c+m-j-1]=comp(tmp,0);
    }
    FFT(B,k,1);
    for(i=0;i<k;i++)B[i]=A[i]*B[i];
    FFT(B,k,-1);
    ans=-1;
    for(i=n-1;i<r;i++)for(j=m-1;j<c;j++){
      tmp=s[i][j]+val;
      if(i>=n)tmp-=s[i-n][j];
      if(j>=m)tmp-=s[i][j-m];
      if(i>=n&&j>=m)tmp+=s[i-n][j-m];
      tmp-=(int(B[i*c+j].r+0.5))*2;
      if(tmp>ans)ans=tmp,sx=i-n+1,sy=j-m+1;
    }
    for(c0=c1=i=0;i<n;i++)for(j=0;j<m;j++)if(a[sx+i][sy+j]==b[i][j])if(b[i][j]=='G')c0++;else c1++;
    printf("Case #%d: %d %d %d %d
",T,sx+1,sy+1,c0,c1);
  }
  return 0;
}