【HDU 5885】XM Reserves

1.题目链接。题目大意：一个n*md的矩阵，每个格点有一个p[i[j],让你选取一个位置,对于所有在以这个点为中心,半径为r的圆内的点的Σp(i,j)/(d+1)最大, d为点到圆心的距离.

2.我们可以对于每个点统计它能对那些格点造成的贡献,然后我们遍历每一个格点就能找到答案了我们从坐标变换的角度来思考,假如一个点的坐标为(x,y),圆心跟它的坐标差为(dx,dy)那么圆心的坐标就是(x+dx,y+dy),由于格点横纵坐标都是整数,我们可以在整数上离散化dx,dy,实际上-r<=dx,dy<=r那么我们对于一个格点,如果它当作圆心(也就是我们选取的位置),剩下能对它产生贡献的点(称为贡献,点)都有一个共同的特性那就是对于每一个贡献点经过一个(dx,dy)的向量偏移后都会到达圆心,即对于所有贡献点(xi+dx,yi+dy)都相等联想到FFT是来求什么的?两个多项式做乘积,能得出结果中每个幂次的系数,我们把每个圆心的坐标看成是多项式乘积结果的每个幂次.其实就是求A(x)=p[i,j],B(x)=1/(sqrt(i*i+j*j)+1.这两个多项式的卷积。

#include <bits/stdc++.h>
using namespace std;
const int maxn = 1 << 21;
const double pi = acos(-1.0);
#define fft FFT
#define r real
struct Complex
{
	double r, i;
	Complex(double _r, double _i) :r(_r), i(_i) {}
	Complex() {}
	Complex operator +(const Complex &b)
	{
		return Complex(r + b.r, i + b.i);
	}
	Complex operator -(const Complex &b)
	{
		return Complex(r - b.r, i - b.i);
	}
	Complex operator *(const Complex &b)
	{
		return Complex(r*b.r - i * b.i, r*b.i + i * b.r);
	}
};
void change(Complex y[], int len)
{
	int i, j, k;
	for (i = 1, j = len / 2; i < len - 1; i++)
	{
		if (i < j)swap(y[i], y[j]);
		k = len / 2;
		while (j >= k)
		{
			j -= k;
			k /= 2;
		}
		if (j < k)j += k;
	}
}
void fft(Complex y[], int len, int on)
{
	change(y, len);
	for (int h = 2; h <= len; h <<= 1)
	{
		Complex wn(cos(-on * 2 * pi / h), sin(-on * 2 * pi / h));
		for (int j = 0; j < len; j += h)
		{
			Complex w(1, 0);
			for (int k = j; k < j + h / 2; k++)
			{
				Complex u = y[k];
				Complex t = w * y[k + h / 2];
				y[k] = u + t;
				y[k + h / 2] = u - t;
				w = w * wn;
			}
		}
	}
	if (on == -1)
		for (int i = 0; i < len; i++)
			y[i].r /= len;
}
Complex a[maxn], b[maxn];
int n, m;
double rr;
int main()
{
	while (~scanf("%d%d%lf", &n, &m, &rr)) {
		int R = ceil(rr);
		int M = max(n, m) + 2 * R;
		int len = 1;
		while (len <= M * M) len <<= 1;
		for (int i = 0; i < len; ++i)
			a[i] = Complex(0.0, 0.0), b[i] = Complex(0.0, 0.0);
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < m; ++j) {
				double p;
				scanf("%lf", &p);
				a[i*M + j] = Complex(p, 0);
			}
		}
		for (int i = -R; i <= R; ++i) {
			for (int j = -R; j <= R; ++j) {
				if (sqrt(i*i + j * j) < rr)
					b[(i + R)*M + j + R] = Complex(1.0 / (sqrt(i*i + j * j) + 1), 0.0);
			}
		}
		FFT(a, len, 1);
		FFT(b, len, 1);
		for (int i = 0; i < len; ++i)
			a[i] = a[i] * b[i];
		FFT(a, len, -1);
		double ans = 0;
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < m; ++j)
				ans = max(ans, a[(i + R)*M + j + R].real);
		}
		printf("%.3lf\n", ans);
	}
	return 0;
}