BZOJ 3527 ZJOI 2014 力 FFT

数轴上有n个点电荷，求各电荷所在位置的电场强度。规定向左为正方向。(强行写题意)

电场力表达式

$$F_j=\sum_{i<j} \frac{q_iq_j}{(i-j)^2}- \sum_{i>j} \frac{q_iq_j}{(i-j)^2}$$
电场强度表达式
$$E_j=\frac{F_j}{q_i}$$
合并可得
$$E_j=\sum_{i<j} \frac{q_i}{(i-j)^2}- \sum_{i>j} \frac{q_i}{(i-j)^2}$$
令 $r_i=\frac{1}{i^2}$，显然有 $r_i=r_{-i}$，且
$$E_j=\sum_{i=0}^{j-1} q_ir_{j-i}-\sum_{i=j+1}^{n-1} q_ir_{i-j}$$
左边就是卷积。右边看起来像卷积。
如果把 $q_i$和 $r_i$绘成一张表，那么左边就是左上-右下对角线的上方，右边就是下方，而下方转过来就是我们想要的了。

Input
5
4006373.885184 15375036.435759 1717456.469144 8514941.004912 1410681.345880
Output
-16838672.693 3439.793 7509018.566 4595686.886 10903040.872

#include <cstdio>
#include <complex>
#include <cmath>
#define rep(i,j,k) for(i=j;i<k;++i)
using namespace std;
typedef complex<double> cd;
const double PI = acos(-1.);
const int N = 500010;
cd q[N], s[N], r[N];
int rev[N];
void fft(cd a[], int n, int ratio) {
    int m, j, k;
    rep(j,0,n) if (j < rev[j]) swap(a[j], a[rev[j]]);
    for (m = 2; m <= n; m *= 2) {
        cd wm(cos(2 * PI / m), sin(ratio * 2 * PI / m));
        for (k = 0; k < n; k += m) {
            cd w(1, 0);
            rep(j, k, k + m / 2) {
                cd u = a[j], t = w * a[j + m / 2];
                a[j] = u + t;
                a[j + m / 2] = u - t;
                w *= wm;
            }
        }
    }
    if (ratio == -1) rep(j,0,n) a[j] /= n;
}
int main() {
    int n, m, i, k;
    double x;
    scanf("%d", &m);
    rep(i,0,m) {
        scanf("%lf", &x);
        q[i] = x; s[m - i - 1] = x;
    }
    rep(i,1,m) r[i] = 1.0 / i / i;
    n = 1; k = -1;
    while (n < 2 * m - 1) n *= 2, ++k;
    rep(i,0,n) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << k);
    fft(q, n, 1); fft(s, n, 1); fft(r, n, 1);
    rep(i,0,n) q[i] *= r[i], s[i] *= r[i];
    fft(q, n, -1); fft(s, n, -1);
    rep(i,0,m) printf("%.3lf\n", q[i].real() - s[m - i - 1].real());
    return 0;
}