HDU5885 XM Reserves（FFT建模）

//-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

HDU - 5885

参考题解：http://blog.csdn.net/crzbulabula/article/details/54577088

题意：给出一个n×m的格子， 格子(i,j)上有值p(i,j)。两个格子之间的距离定义为两个格子中心的欧几里得距离（sqrt(dx*dx+dy*dy)）。你可以选择一个格子， 距离它小于等于r的那些格子都会贡献p(i,j)/(d+1)，d为两个格子的距离。选择一个格子，使得贡献和最大。

题解：

两个点之间的欧几里得距离是sqrt(dx*dx+dy*dy)，必定与四个量（x1,y1,x2,y2）有关。

如果要用FFT求解，必须简化四个量。因此我们枚举满足条件的偏移向量(dx,dy)。

偏移向量要怎么用呢？我们可以把(x,y)+(dx,dy)=(x+dx,y+dy)这个性质应用到指数上面去。

(x,y)的价值是p(x,y)，它经过偏移(dx,dy)会对(x+dx,y+dy)贡献p(i,j)/( sqrt(dx*dx+dy*dy) +1)。

现在问题已经简化了，但是它是两维的，我们要用FFT必须把它简化成一维。

A[i∗M+j]=p(i,j)，B[dx * M + dy] = 1/( sqrt(dx*dx+dy*dy) +1)​​ ，它们的卷积 C = A * B 的第 (i * M + j) 项就是 (i, j) 格子的贡献和。

但是dx,dy可能是负数，所以我们把所有的dx替换成dx+R

最后的求解过程如下：

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;

const double PI=acos(-1);
struct C {
	double x,y;
	C(double _x=0,double _y=0) {
		x=_x;y=_y;
	}
	C operator+(const C &p) {
		return C(x+p.x,y+p.y);
	}
	C operator-(const C &p) {
		return C(x-p.x,y-p.y);
	}
	C operator*(const C &p) {
		return C(x*p.x-y*p.y,x*p.y+y*p.x);
	}
};
void fft(C x[],int len,int on) {
	for(int i=1,j=len/2;j<len-1;++i) {
		if(i<j) swap(x[i],x[j]);
		int k=len>>1;
		while(j>=k) {
			j-=k;
			k>>=1;
		}
		if(j<k) j+=k;
	}
	for(int h=2;h<=len;h<<=1) {
		C wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
		for(int j=0;j<len;j+=h) {
			C w(1,0);
			for(int k=j;k<j+h/2;++k) {
				C u=x[k];
				C t=w*x[k+h/2];
				x[k]=u+t;
				x[k+h/2]=u-t;
				w=w*wn;
			}
		}
	}
	if(on==-1) {
		for(int i=0;i<len;++i) x[i].x/=len;
	}
}

const int N=(1<<21);
C X1[N],X2[N];

double dis(int x,int y) {
    return sqrt(x*x+y*y);
}

int main() {
    int n,m;double r;
    while(~scanf("%d%d%lf",&n,&m,&r)) {
        
        int R=ceil(r);
        int M=max(n+2*R,m+2*R);
        
        int len=1;
        while(len<=M*M) len<<=1;
        
        for(int i=0;i<len;++i) X1[i]=X2[i]=C(0,0);
        for(int i=0;i<n;++i) {
            for(int j=0;j<m;++j) {
                int x;scanf("%d",&x);
                X1[i*M+j]=C(x,0);
            }
        }
        for(int i=-R;i<=R;++i) {
            for(int j=-R;j<=R;++j) {
                if(dis(i,j)<r) X2[(i+R)*M+j+R]=C(1.0/(dis(i,j)+1),0);
            }
        }
        fft(X1,len,1);
        fft(X2,len,1);
        for(int i=0;i<len;++i) X1[i]=X1[i]*X2[i];
        fft(X1,len,-1);
        double ans=0;
        for(int i=0;i<n;++i) {
            for(int j=0;j<m;++j) {
                ans=max(ans,X1[(i+R)*M+j+R].x);
            }
        }
        printf("%.3lf\n",ans);
    }
    return 0;
}