XM Reserves(FFT建模)

描述

传送门:2016 ACM/ICPC Asia Regional Qingdao Online- Problem H

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$ ($1≤n,m≤500$) and one float-point number $r$ ($0≤r≤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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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

1
2
3
9.000
24.142
17.956

思路

  • 直接暴力算的话复杂度是$O(nmr^2)$显然是会超时的。
  • 若利用FFT计算,可以用一个向量表示位置信息和值,另一个向量表示某一个位置平移后的权。这样的复杂度为$O(max(mn,r ^ 2)\log(\max(mn,r ^ 2)))$。
  • 向量A表示位置信息的值,最大为$m \times n$;向量B表示某一个偏移量下的权,最大为$r\times r$。向量C表示卷积得到的向量。
    $$C(i,j) = \sum_{i_1+i_2=i}\sum_{j_1+j_2=j} A(i_1,j_1) \times B(i_2,j_2)$$
    注意到向量C最大为$((m+2R) \times (n+2R))$,方便起见,将其转化为$(M \times M)$的向量,其中$M=\max(m,n)+2R$。
  • 这样就可以把二维的情况映射到一维了。
    $$A(i,j) = A[i \times M + j]$$
    $$B(i,j) = B[(i+R) \times M + j + R]$$
  • 这里加上$R$是因为$B$中的偏移量是在$[-R, R]$,但是实现的时候我们不能使用负下标。
  • 最后的结果就是向量C:$$C(i,j) = C[M \times (i+R)+j+R]$$

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <bits/stdc++.h>
using namespace std;

const double PI = acos(-1.0);

struct Complex
{
double x, y;
Complex(double x = 0.0, double y = 0.0) : x(x), y(y) {}
Complex operator-(const Complex& b) const { return Complex(x - b.x, y - b.y); }
Complex operator+(const Complex& b) const { return Complex(x + b.x, y + b.y); }
Complex operator*(const Complex& b) const { return Complex(x * b.x - y * b.y, x * b.y + y * b.x); }
};

void change(Complex y[], int len)
{
for (int i = 1, j = len / 2; i < len - 1; i++)
{
if (i < j) swap(y[i], y[j]);
int 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].x /= len;
}

const int N = 1 << 20;
Complex a[N], b[N];

int main()
{
int n, m;
double r;
while (~scanf("%d%d%lf", &n, &m, &r))
{
memset(a, 0, sizeof(a));
memset(b, 0, sizeof(b));
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 < n; i++)
for (int j = 0; j < m; j++)
scanf("%lf", &a[i * M + j].x), a[i * M + j].y = 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] = {1.0 / (sqrt(i * i + j * j) + 1), 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].x);
printf("%.3f\n", ans);
}
}
捐助作者
0%