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\times b$即可。

时间复杂度$O(BRC\log(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\n",T,sx+1,sy+1,c0,c1);

  }

  return 0;

}