HDU5885 XM Reserves (FFT)

HDU-5885 (FFT)

可以看到题目是一个二维的作差转移，同理的，我们可以将二维转移转化为序列\((x,y)\rightarrow x\cdot m+y\)，但是要注意转移边界的问题，建议在每一行多加一些，即\((x,y)\rightarrow x\cdot 3m+y+m\)

如果你还不会作差卷积的构造: BZOJ-3527 题解

处理完转移之后用\(FFT\)即可

#include<bits/stdc++.h>
using namespace std;

//#define double long double

#define reg register
typedef long long ll;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)

template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); } 
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); } 

char IO;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const double PI=acos(-1);
const int N=(1<<21)+5;

int n,m;
double r;
struct Cp{
	double x,y;
	Cp(){}
	Cp(double _x,double _y){ x=_x,y=_y; }
	Cp operator + (const Cp t){ return Cp(x+t.x,y+t.y); }
	Cp operator - (const Cp t){ return Cp(x-t.x,y-t.y); }
	Cp operator * (const Cp t){ return Cp(x*t.x-y*t.y,x*t.y+y*t.x); }
};
Cp a[N],b[N];
int rev[N];

void FFT(int n,Cp *a,int f){
	rep(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(reg int i=1;i<n;i<<=1) {
		Cp w(cos(PI/i),f*sin(PI/i));
		for(reg int l=0;l<n;l+=i*2) {
			Cp e(1,0);
			for(reg int j=l;j<l+i;j++,e=e*w) {
				Cp t=a[j+i]*e;
				a[j+i]=a[j]-t;
				a[j]=a[j]+t;
			}
		}
	}
	if(f==-1) rep(i,0,n-1) a[i].x/=n;
}


int main(){
	while(scanf("%d%d%lf",&n,&m,&r)==3) {
		rep(i,0,n-1) rep(j,0,m-1) scanf("%lf",&a[i*(m*3)+j+m].x);//初始,偏移为m*3
        //转移
		rep(i,-n+1,n-1) rep(j,-m+1,m-1) if(i*i+j*j<r*r) b[i*m*3+j+3*n*m].x+=1.0/(sqrt(i*i+j*j)+1);
		int R=1,c=0;
		while(R<=n*m*6) R<<=1,c++;//防止转移越界，开到6倍
		rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<(c-1));
		FFT(R,a,1),FFT(R,b,1);
		rep(i,0,R) a[i]=a[i]*b[i];
		FFT(R,a,-1);
		double ans=0;
		rep(i,0,R) {
			int x=(i-3*n*m)/(3*m),y=(i-3*n*m)%(3*m)-m;
			if(x>=0 && x<n && y>=0 && y<m) cmax(ans,a[i].x);
		} // 计算答案
		rep(i,0,R) b[i].x=b[i].y=a[i].x=a[i].y=0;
		printf("%.3lf\n",ans);
	}
}