「Bzoj 3527」「ZJOI2014」力 (卷积FFT)

BZOJ 3527
给定$n, q_i$, 求
$$
E_j=\sum\limits_{i=1}^{j-1} \frac{q_i}{(i - j)^2} - \sum\limits_{i=j+1}^{n} \frac{q_i}{(i - j)^2}
$$

设$f(i)=q_i, g(i)=\frac 1{i^2}$
则
$$
E_j=\sum\limits_{i=1}^{j-1} f(i)g(i-j) - \sum\limits_{i=j+1}^{n} f(i)g(i-j)
$$
因为$g(i)$为偶函数，则
$$
E_j=\sum\limits_{i=1}^{j-1} f(i)g(j-i) - \sum\limits_{i=j+1}^{n} f(i)g(i-j)
$$
那么前面是一个卷积形式，后面可以翻转$f(i)$后变为另一个卷积形式，FFT 计算即可。

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<map>
#include<vector>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky {

    const int MAXN = 300000 + 5;
    const db PI = acos(-1);

    struct C {
        db r, i;
        C() {r = i = 0;}
        C(db x, db y) {r = x, i = y;}
        C operator + (C &x) {return C(r + x.r, i + x.i);}
        C operator - (C &x) {return C(r - x.r, i - x.i);}
        C operator * (C &x) {return C(r * x.r - i * x.i, x.i * r + x.r * i);}
        void operator += (C &x) {r += x.r, i += x.i;}
        void operator *= (C &x) {
            db tr = r, ti = i;
            r = tr * x.r - ti * x.i;
            i = x.i * tr + x.r * ti;
        }
    } f[MAXN], g[MAXN];

    int m, n, r[MAXN];
    db p[MAXN], q[MAXN], a[MAXN], b[MAXN], gg[MAXN];

    void FFT(C *a, int op) {
        for (int i = 0; i < n; ++i)
            if (i < r[i]) swap(a[i], a[r[i]]);
        for (int i = 1; i < n; i <<= 1) {
            C Wn = C(cos(PI / i), sin(PI / i) * op);
            for (int j = 0; j < n; j += (i << 1)) {
                C w = C(1, 0), *a0 = a + j, *a1 = a0 + i;
                for (int k = 0; k < i; ++k) {
                    C t = *a1 * w;
                    *a1 = *a0 - t, *a0 += t;
                    ++a0, ++a1, w *= Wn;
                }
            }
        }
    }
    void pro(db *tf, db *tg, db *ret) {
        for (int i = 0; i < n; ++i) 
            f[i].r = tf[i], g[i].r = tg[i], f[i].i = 0, g[i].i = 0;
        FFT(f, 1), FFT(g, 1);
        for (int i = 0; i < n; ++i) f[i] *= g[i];
        FFT(f, -1);
        for (int i = 0; i <= m; ++i) ret[i] = (fabs(f[i].r) / n);
    }

    void clean() {
        ms(r, 0);
    }
    int solve() {

        clean();
        scanf("%d", &n);
        for (int i = 1; i <= n; ++i) {
            scanf("%lf", &p[i]);
            q[i] = p[i], gg[i] = 1.0 / i / i;
        }

        reverse(q + 1, q + 1 + n);

        int l = 0;
        for (m = n, n = 1; n <= m + m; n <<= 1, ++l);

        for (int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));

        pro(p, gg, a);
        pro(q, gg, b);

        for (int i = 1; i <= m; ++i) printf("%.5f\n", a[i] - b[m - i + 1]);

        return 0;
    }
}
int main() {
    flyinthesky::solve();
    return 0;
}