题目地址：洛谷：【P3338】[ZJOI2014]力 – 洛谷、BZOJ：Problem 3527. — [Zjoi2014]力

题目描述

给出n个数q_i，给出F_j的定义如下：
F_j = \sum_{i < j} \frac{q_iq_j}{(i-j)^2} - \sum_{i > j} \frac{q_iq_j}{(i-j)^2}
令E_i = \frac{F_i}{q_i}，求E_i。

输入输出格式

输入格式：
第一行一个整数n。
接下来n行每行输入一个数，第i行表示qi。

输出格式：
n行，第i行输出Ei。
与标准答案误差不超过1e-2即可。

输入输出样例

输入样例#1：

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

输出样例#1：

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

说明

对于30%的数据，n≤1000。
对于50%的数据，n≤60000。
对于100%的数据，n≤100000，0<qi<1000000000。

题解

根据F的定义我们展开E的定义，如下
E_i = \sum_{j=1}^{i-1} \frac{q_j}{(i-j)^2} - \sum_{j=i+1}^{n} \frac{q_j}{(j-i)^2}
如果我们让a_i = q_i, b_i = \frac{1}{i^2}, b_0 = 0，代换掉再观察E的定义
E_i = \sum_{j=0}^{i-1} a_jb_{i-j} - \sum_{j=i+1}^{n} a_jb_{j-i}
第一项的求和实际上已经是卷积的形式了，但是后面这一项不是。我们考虑给这一项改一改
E_i = \sum_{j=0}^{i-1} a_jb_{i-j} - \sum_{j=0}^{n-i-1} a_{i+j}b_{j}
后面这一项好像还不是卷积呀，那如果我们令c_{n-i-j-1}=a_{i+j}（相当于把a数组翻过来存），带换掉？
E_i = \sum_{j=0}^{i-1} a_jb_{i-j} - \sum_{j=0}^{n-i-1} c_{n-i-j-1}b_{j}
好像可以卷积了！为了便于理解，我们换一种形式
A_i = \sum_{j=0}^{n-i-1} c_{n-i-j-1}b_{j} \Rightarrow A_{n-i-1} = \sum_{j=0}^{i} c_{i-j}b_{j}
这样我们看后面一项就是个卷积的形式。

代码

// Code by KSkun, 2018/3
#include <cstdio>
#include <cmath>
#include <cstring>

#include <algorithm>

const int MAXN = 1 << 20;
const double PI = std::acos(-1);

struct Complex {
    double real, imag;
    Complex(double real = 0, double imag = 0) : real(real), imag(imag) {}
    inline Complex operator+(const Complex &rhs) const {
        return Complex(real + rhs.real, imag + rhs.imag);
    }
    inline Complex operator-(const Complex &rhs) const {
        return Complex(real - rhs.real, imag - rhs.imag);
    }
    inline Complex operator*(const Complex &rhs) const {
        return Complex(real * rhs.real - imag * rhs.imag, real * rhs.imag + imag * rhs.real);
    }
    inline Complex& operator*=(const Complex &rhs) {
        return *this = *this * rhs;
    }
};

int n, m, len, rev[MAXN];
Complex a[MAXN], b[MAXN];
double s[MAXN], ans[MAXN];

inline void fft(Complex *arr, int f) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) std::swap(arr[i], arr[rev[i]]);
    }
    for(int i = 1; i < n; i <<= 1) {
        Complex wn(std::cos(PI / i), f * std::sin(PI / i));
        for(int j = 0; j < n; j += i << 1) {
            Complex w(1, 0);
            for(int k = 0; k < i; k++) {
                Complex x = arr[j + k], y = w * arr[j + k + i];
                arr[j + k] = x + y;
                arr[j + k + i] = x - y;
                w *= wn;
            }
        }
    }
}

int N;

int main() {
    scanf("%d", &N); n = N - 1;
    for(int i = 0; i < N; i++) {
        scanf("%lf", &s[i]);
    }
    m = n << 1;
    for(n = 1; n <= m; n <<= 1) len++;
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
    }
    for(int i = 0; i < N; i++) {
        a[i].real = s[i];
        if(i) b[i].real = 1 / double(i) / double(i);
    }
    fft(a, 1);
    fft(b, 1);
    for(int i = 0; i <= n; i++) {
        a[i] *= b[i];
    }
    fft(a, -1);
    for(int i = 0; i < N; i++) {
        ans[i] = a[i].real / double(n);
    }
    memset(a, 0, sizeof(a));
    for(int i = 0; i < N; i++) {
        a[i].real = s[N - i - 1];
    }
    fft(a, 1);
    for(int i = 0; i <= n; i++) {
        a[i] *= b[i];
    }
    fft(a, -1);
    for(int i = 0; i < N; i++) {
        ans[i] -= a[N - i - 1].real / double(n);
    }
    for(int i = 0; i < N; i++) {
        printf("%.3lf\n", ans[i]);
    }
    return 0;
}