hdu 5885 XM Reserves（FFT）

Posted 2023-02-22 jpphy0

目录

• 问题
• 分析
• 代码

问题

hdu 5885 XM Reserves - https://acm.hdu.edu.cn/showproblem.php?pid=5885

分析

• 二维数组地址线性化
• 一个格点对周边格点的产生影响，被影响的格点的地址从线性地址角度看，可以用多项式乘积的幂来表达
• 对值的影响，可以用多项式的系数来表达
• 格点$(i,j)$对格点$(i+R-x,j+R-y)$的贡献：
  $$\Delta v=v_{i,j}\cdot z^{i\cdot (M+2R)+j}\times \frac{1}{d+1}\cdot z^{(R-x)\cdot (M+2R)+R-y}，当d=\sqrt{x^2+y^2}<R$$
  

代码

#include<bits/stdc++.h>
using namespace std;
typedef complex<double> cd;
const double DFT = 2.0, IDFT  = -2.0, PI = acos(-1);
const int MX = 1100+10, MXL = MX*MX;
int n, m, rev[MXL];
cd a[MXL], b[MXL];
void fft(cd p[], int lim, int r[], double mode)
    for(int i = 0; i < lim; ++i) if(i < r[i]) swap(p[i], p[r[i]]);
    for(int len = 2; len <= lim; len <<= 1)
        cd wn(cos(mode*PI/len), sin(mode*PI/len));
        for(int s = 0; s < lim; s += len)
            cd w(1.0, 0.0);
            for(int cur = s; cur < s+(len>>1); ++cur, w *= wn)
                cd l = p[cur], r = w*p[cur+(len>>1)];
                p[cur] = l + r, p[cur+(len>>1)] = l - r;
            
        
    
    if(mode == DFT) return;
    for(int i = 0; i < lim; ++i) p[i] /= lim;

int main()
    double r, tmp, ans;    
    while(scanf("%d%d%lf", &n, &m, &r) == 3)
        int R = (int)(r+0.5), M = m+(R<<1);
        int lim = M*(n + (R << 1)), bit = 1;
        while(lim >>= 1) ++bit;
        lim = 1 << bit;
        for(int i = 0; i < lim; ++i) rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
        for(int i = 0; i < lim; ++i) a[i] = b[i] = 0;
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < m; ++j) scanf("%lf", &tmp), a[i*M+j] = tmp;
        for(int i = 0; i <= R<<1; ++i)
            for(int j = 0; j <= R<<1; ++j)
                tmp = sqrt((R-i)*(R-i) + (R-j)*(R-j));
                if(tmp < r ) b[i*M+j] = 1.0/(tmp+1.0);
                    
        fft(a, lim, rev, DFT), fft(b, lim, rev, DFT);
        for(int i = 0; i < lim; ++i) a[i] *= b[i];
        fft(a, lim, rev, IDFT);
        ans = 0;
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < m; ++j) ans = max(a[(i+R)*M+j+R].real(), ans);
        printf("%.3lf\\n", ans);
    
    return 0;