BZOJ 3527. [Zjoi2014]力

黑暗爆炸上没有spj。

第一部分的答案为 $\sum \limits_{j = 1}^{i - 1} q_j \times \frac{1}{(i-j)^2}$，直接卷积就行。
第二部分就把 q 数组reverse一下再卷一遍。

#include <bits/stdc++.h>
#define db long double

const int N = 100000 * 3 + 7;
db a[N], ans[N][2];

namespace FFT {
    const db pi = acos(-1.0);
    struct Complex {
        db r, i;
        Complex() {}
        Complex(db r, db i): r(r), i(i) {}
        void clear() { r = i = 0.0; }
        Complex operator + (const Complex &p) const { return Complex(r + p.r, i + p.i); }
        Complex operator - (const Complex &p) const { return Complex(r - p.r, i - p.i); }
        Complex operator * (const Complex &p) const { return Complex(r * p.r - i * p.i, r * p.i + i * p.r); }
    } A[N], B[N];
    int n, l, r[N];
    void init(int m) {
        n = 1, l = 0;
        while (n <= m) n <<= 1, l++;
        for (int i = 0; i < n; i++) {
            r[i] = r[i >> 1] >> 1 | ((i & 1) << (l - 1));
        }
    }
    void FFT(Complex *a, int pd) {
        for (int i = 0; i < n; i++)
            if (i < r[i])
                std::swap(a[i], a[r[i]]);
        for (int mid = 1; mid < n; mid <<= 1) {
            Complex wn(cos(pi / mid), pd * sin(pi / mid));
            for (int l = mid << 1, j = 0; j < n; j += l) {
                Complex w(1.0, 0.0);
                for (int k = 0; k < mid; k++, w = w * wn) {
                    Complex u = a[k + j], v = w * a[k + j + mid];
                    a[k + j] = u + v;
                    a[k + j + mid] = u - v;
                }
            }
        }
        if (pd == -1)
            for (int i = 0; i < n; i++)
                a[i] = Complex(a[i].r / n, a[i].i / n);
    }
    void Mul(Complex *a, Complex *b) {
        FFT(a, 1); FFT(b, 1);
        for (int i = 0; i < n; i++)
            a[i] = a[i] * b[i];
        FFT(a, -1);
    }
    void solve(const db *a, int opt, int m) {
        for (int i = 0; i < n; i++)
            A[i].clear(), B[i].clear();
        for (int i = 1; i <= m; i++)
            A[i] = Complex(a[i], 0.0);
        for (int i = 1; i <= m; i++)
            B[i] = Complex(1.0 / (1.0 * i) / (1.0 * i), 0.0);
        Mul(A, B);
        for (int i = 1; i <= m; i++)
            ans[i][opt] = A[i].r;
    }
}

int n;

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)
        scanf("%Lf", a + i);
    FFT::init(n * 2);
    FFT::solve(a, 1, n);
    std::reverse(a + 1, a + 1 + n);
    FFT::solve(a, 0, n);
    for (int i = 1; i <= n; i++)
        printf("%.3Lf\n", ans[i][1] - ans[n - i + 1][0]);
    return 0;
}