HDU5885 XM Reserves（FFT建模）

来源：互联网发布：怎么在淘宝里卖二东西编辑：程序博客网时间：2024/05/21 09:47

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

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;}

0 0