ZJOI 2014 力 FFT 快速傅里叶变换

https://www.luogu.org/problemnew/show/P3338
这是学了$FFT$后除了板子的一题QAQ
这种题目没做过不知道怎么下手看了题解才知道

记$i < j$ 时 $f(x) = \frac{1}{x^2}$, $i > j$时 $f(x) = \frac{-1}{x^2}$否则$f(x) = 0$
则$E_i = \sum_{j = 0}^{n - 1}f(i - j)×q_j$
我们发现$f$函数下标可能是负的，平移$n$个单位即可

对于卷积直接$FFT$优化就好啦
题目做少了以后还是多做一点把

#include<bits/stdc++.h>

using namespace std;

const int N = 1e6 + 10;
const double Pi = acos(-1);

struct Complex {
    double x, y;
    Complex(double xx = 0, double yy = 0) {
        x = xx, y = yy;
    }
    Complex operator + (const Complex &T) {
        return Complex(x + T.x, y + T.y);
    }
    Complex operator - (const Complex &T) {
        return Complex(x - T.x, y - T.y);
    }
    Complex operator * (const Complex &T) {
        return Complex(x * T.x - y * T.y, x * T.y + y * T.x);
    }
}a[N], b[N];

void FFT_Init(int limit, Complex *a) {
    int k = log2(limit);
    for(int i = 0; i < limit; ++ i) {
        int t = 0;
        for(int j = 0; j < k; ++ j)
            if((1 << j) & i)
                t |= 1 << (k - j - 1);
        if(t < i) swap(a[t], a[i]);
    }
}

void FFT(int n, Complex *a, int fh) {
    FFT_Init(n, a);
    for(int limit = 2; limit <= n; limit <<= 1) {
        double xita = 2.0 * Pi / limit;
        Complex Wn = Complex(cos(xita), sin(xita) * fh), W = Complex(1, 0);
        for(int j = 0; j < n; j += limit, W = Complex(1, 0)) 
            for(int i = j; i < j + (limit >> 1); ++ i, W = W * Wn) {
                Complex a1 = a[i], a2 = a[i + (limit >> 1)] * W;
                a[i] = a1 + a2, a[i + (limit >> 1)] = a1 - a2;
            }
    }
    if(fh == -1)
        for(int i = 0; i < n; ++ i)
            a[i].x /= n;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("3338.in", "r", stdin);
    freopen("3338.out", "w", stdout);
#endif
    int n, limit = 1;
    scanf("%d", &n);
    while(limit <= (n - 1) * 3)
        limit <<= 1;
    for(int i = 0; i < n; ++ i)
        scanf("%lf", &b[i].x);
    for(int i = 0; i < (n << 1) - 1; ++ i) {
        a[i].x = i - n + 1;
        if(a[i].x != 0) 
            a[i].x = (i < n ? -1.0 : 1.0) / (a[i].x * a[i].x);
    }
    FFT(limit, a, 1), FFT(limit, b, 1);
    for(int i = 0; i < limit; ++ i)
        a[i] = a[i] * b[i];
    FFT(limit, a, -1);
    for(int i = n - 1; i < (n << 1) - 1; ++ i)
        printf("%.6lf\n", a[i].x);
    return 0;
}
// 歌声还在游走 你榴花般的双眸