题意:
给出一个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;
}