HDU-5885 XM Reserves

题面

Description

As an eligible Ingress Resistance Agent you should know your power source, the Exotic Matter.
We call it XM, which is the driving force behind all of our actions in Ingress.
XM allows us to construct items through hacking portals, to attack enemy portals, make links and create fields.

We try to collect XM from the ground. XM concentration come from location based services, meaning that areas with a lot of foot traffic have higher amounts versus places that don't.
You can collect XM by moving through those areas.
The XM will be automatically harvested by your Scanner when it is within your interaction circle/range.

Alice decides to select a location such that she can collect XM as much as possible.
To simplify the problem, we consider the city as a grid map with size `n*m' numbered from ((0,0)) to ((n-1,m-1)).
The XM concentration inside the block ((i,j)) is (p(i,j)).
The radius of your interaction circle is (r).
We can assume that XM of the block ((i,j)) are located in the centre of this block.
The distance between two blocks is the Euclidean distance between their centres.

Alice stands in the centre of one block and collects the XM.
For each block with the distance d smaller than r to Alice, and whose XM concertation is (p(i,j)), Alice's scanner can collects (p(i,j)/(1+d)) XM from it.

Help Alice to determine the maximum XM which she can collect once he stands in the centre of one block.

Input

There are multiple cases.
For each case, the first line consists two integers (n, m~(1le n,m le 500)) and one float-point number (r~(0le r le 300)).
Each of the following (n) line consists (m) non-negative float-point numbers corresponding to the XM concentrations inside each blocks.

Output

For each case, output the maximum XM which Alice can collect in one line.
Your answers should be rounded to three decimal places.

Sample Input

3 3 1
1 3 6
7 9 4
2 8 1
3 3 2
1 3 6
7 9 4
2 8 1
5 5 1.5
4 3 2 9 1 
3 4 3 2 8
9 4 3 2 1
2 3 0 1 2
6 3 4 3 1

Sample Output

9.000
24.142
17.956

题意

给定一个方格图，每格有一定数量的道具(p(i,j))，有一个工具，可以吸取半径为r范围内的道具，假设人物所处位置是方格的中心，所有的道具也位于方格的中心，两点之间的距离为欧式距离。当所在点与某格之间的距离d小于r时，可以吸收(p(i,j)/(1+d))数量的道具，选定一个初始点，使吸收的道具数量最大。保留3位小数

题解

直接扫一遍，显然复杂度爆炸，但我们想要获取答案，就要知道每一格能吸收多少的道具。怎么避免扫一遍呢？

我们可以考虑偏移量。

对于一个点，能对其产生贡献的点，是由很多其他的点((x,y))组成的，而这些点加上一个偏移量((dx,dy))，就是该点，这个偏移量显然要满足((d =sqrt{dx^2+dy^2}) < r)，其产生的贡献就是(p(i,j)/(1+d))，这个形式，是不是有点fft的意思了

我们把每个点的道具数量(p(i,j))放进一个数组a，对应的偏移量的贡献(1/(1+d))放进另一个数组b，进行卷积，就得到了每个点的答案。也就是不同点，进行不同大小的偏移，到达了某个点。在fft形式中，也就是每个点对应了一个指数，其系数就是答案。但是由于fft是一维的，我们的偏移量和p都是二维的，所以我们要把二维化成一维，这个很简单，直接用(i*M+j)表示就可以了，(M=max(n,m)+2*R)，另外，偏移量有正有负，所以我们要把行偏移量和列偏移量都加上一个R，防止负数出现，这样卷积后，对应的(a[(i+R)*M+j+R].r)就是答案。i和j都加了R是因为b的((dx,dy))都加了R。

这样我们扫一遍系数，取最大值即可

代码

#include <bits/stdc++.h>
using namespace std;
const double pi = acos(-1.0);
const int N = 2e6 + 50;
typedef long long ll;
struct cp {
    double r, i;
    cp(double r = 0, double i = 0): r(r), i(i) {}
    cp operator + (const cp &b) {
        return cp(r + b.r, i + b.i);
    }
    cp operator - (const cp &b) {
        return cp(r - b.r, i - b.i);
    }
    cp operator * (const cp &b) {
        return cp(r * b.r - i * b.i, r * b.i + i * b.r);
    }
};
void change(cp a[], int len) {
    for (int i = 1, j = len / 2; i < len - 1; i++) {
        if (i < j) swap(a[i], a[j]);
        int k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) j += k; 
    }
}
void fft(cp a[], int len, int op) {
    change(a, len);
    for (int h = 2; h <= len; h <<= 1) {
        cp wn(cos(-op * 2 * pi / h), sin(-op * 2 * pi / h));
        for (int j = 0; j < len; j += h) {
            cp w(1, 0);
            for (int k = j; k < j + h / 2; k++) {
                cp u = a[k];
                cp t = w * a[k + h / 2];
                a[k] = u + t;
                a[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if (op == -1) {
        for (int i = 0; i < len; i++) a[i].r /= len;
    }
}
cp a[N], b[N];
double p[550][550];
int main() {
    int n, m;
    double r;
    while (~scanf("%d%d%lf", &n, &m, &r)) {
        int R = ceil(r);
        
        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] = cp(0, 0);
            b[i] = cp(0, 0);
        }
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                scanf("%lf", &p[i][j]);
                a[i * M + j] = cp(p[i][j], 0);
            }
        }
        for (int i = -R; i <= R; i++) {
            for (int j = -R; j <= R; j++) {
                if (sqrt(i * i + j * j) < r) {
                    b[(i + R) * M + j + R] = cp(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.0;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                ans = max(ans, a[(i + R) * M + j + R].r);
            }
        }
        printf("%.3f
", ans);
    }
    return 0;
}