洛谷 P3723 [AH2017/HNOI2017]礼物

由题面可得：

\[E_j = \sum_{i = 1}^{j - 1} \frac{q_i}{(i - j)^2} - \sum_{i = j + 1}^{n} \frac{q_i}{(i - j)^2} \]

令 \(q_0 = 0\)，并将没有意义的分式的值视为 \(0\)，则有：

\[E_j = \sum_{i = 0}^j \frac{q_i}{(i - j)^2} - \sum_{i = j}^n \frac{q_i}{(i - j)^2} \]

令 \(A(i) = q_i, B(i) = \dfrac 1{i^2}\)，则原式可化为：

\[E_j = \sum_{i = 0}^j[A(i) \times B(j - i)] - \sum_{i = 0}^{n - j}[A(i + j) \times B(i)] \]

令 \(A'(i) = A(n - i)\)，则有：

\[\begin{aligned} E_j &= \sum_{i = 0}^j[A(i) \times B(j - i)] - \sum_{i = 0}^{n - j}[A'(n - i - j) \times B(i)] \\ &= (A * B)[j] - (A' * B)[n - j] \end{aligned} \]

FFT 加速卷积即可。

代码：

#include <bits/stdc++.h>

#define MAXN 400100

using namespace std;

const double PI = acos(-1);

int n, len = 1, bits, rev[MAXN];
double q[MAXN];

struct Complex {
    double r, i;

    Complex operator+(const Complex &rhs) const {
        return {r + rhs.r, i + rhs.i};
    }

    Complex operator-(const Complex &rhs) const {
        return {r - rhs.r, i - rhs.i};
    }

    Complex operator*(const Complex &rhs) const {
        return {r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r};
    }
} a[MAXN], b[MAXN], a1[MAXN];

void FFT(Complex A[], int flag) {
    for (int i = 0; i < len; i++) if (rev[i] > i) swap(A[i], A[rev[i]]);
    for (int i = 1; i < len; i <<= 1) {
        Complex wn = {cos(PI / i), flag * sin(PI / i)};
        for (int j = 0; j < len; j += (i << 1)) {
            Complex w = {1, 0};
            for (int k = j; k < j + i; k++) {
                Complex t = w * A[i + k];
                A[k + i] = A[k] - t;
                A[k] = A[k] + t;
                w = w * wn;
            }
        }
    }
}

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) scanf("%lf", &q[i]);
    for (int i = 1; i <= n; i++) a[i].r = a1[n - i].r = q[i];
    for (int i = 1; i <= n; i++) b[i].r = 1.0 / i / i;
    while (len < (n << 1)) len <<= 1, bits++;
    for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bits - 1));
    FFT(a, 1), FFT(a1, 1), FFT(b, 1);
    for (int i = 0; i < len; i++) a[i] = a[i] * b[i], a1[i] = a1[i] * b[i];
    FFT(a, -1), FFT(a1, -1);
    for (int i = 1; i <= n; i++) printf("%lf\n", (a[i].r - a1[n - i].r) / len);
    return 0;
}