hdu5885

题意：

给出一个n*m的网格，每个格点有一个权值p(i,j)，每个格点可以对周围距离严格小于r的格点贡献权值除以欧几里得距离+1，问最后所有点中价值最大的是多少

solution：

每个格点可以贡献的范围显然是一个圆，令R = ceil(r)，将原网格扩充为N = n + 2*R，M = m + 2*R，那么，贡献的转移就像是二维空间上的卷积一样，将二维上的点和转移方式转换成一维，FFT解决

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<vector>
#include<queue>
#include<set>
#include<map>
#include<stack>
#include<bitset>
#include<ext/pb_ds/priority_queue.hpp>
using namespace std;

const int maxn = 3E6 + 20;
typedef double DB;
const DB eps = 1E-9;
const DB PI = acos(-1);

struct Virt{
	DB r,i; Virt(){}
	Virt(DB r,DB i): r(r),i(i){}
	Virt operator + (const Virt &b) {return Virt(r + b.r,i + b.i);}
	Virt operator - (const Virt &b) {return Virt(r - b.r,i - b.i);}
	Virt operator * (const Virt &b) {return Virt(r*b.r - i*b.i,r*b.i + i*b.r);}
}A[maxn],B[maxn],C[maxn];

int n,m,R,M,N;
DB r;

void Rader(Virt *F)
{
	int j = (N >> 1);
	for (int i = 1; i < N - 1; i++)
	{
		if (i < j) swap(F[i],F[j]);
		int k = (N >> 1);
		while (j >= k) j -= k,k >>= 1;
		j += k;
	}
}

void FFT(Virt *F,int on)
{
	Rader(F); DB G = 2.00*PI*(DB)(on);
	for (int k = 2; k <= N; k <<= 1)
	{
		Virt wn = Virt(cos(G / (DB)(k)),sin(G / (DB)(k)));
		for (int i = 0; i < N; i += k)
		{
			Virt w = Virt(1,0);
			for (int j = i; j < i + (k >> 1); j++)
			{
				Virt u = F[j],t = w * F[j + (k >> 1)];
				F[j] = u + t; F[j + (k >> 1)] = u - t;
				w = w * wn;
			}
		}
	}
	if (on == -1)
		for (int i = 0; i < N; i++) F[i].r /= (DB)(N);
}

int main()
{
	
	while (scanf("%d%d%lf",&n,&m,&r) != EOF)
	{
		R = ceil(r); M = m + 2*R; N = 1;
		while (N < (n + 2*R)*(m + 2*R)) N <<= 1;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
			{
				DB x; scanf("%lf",&x);
				A[i*M + j] = Virt(x,0);
			}
		for (int i = -R; i <= R; i++)
			for (int j = -R; j <= R; j++)
			{
				DB dis = sqrt(i*i + j*j);
				if (r > dis)
					B[(i + R)*M + j + R] = Virt(1.00 / (dis + 1.00),0);
			}
		FFT(A,1); FFT(B,1);
		for (int i = 0; i < N; i++) C[i] = A[i] * B[i];
		FFT(C,-1); DB Ans = 0;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				Ans = max(Ans,C[(i + R)*M + j + R].r);
		printf("%.3f\n",Ans);
		for (int i = 0; i < N; i++) A[i] = B[i] = C[i] = Virt(0.00,0.00);
	}
	return 0;
}