Link:
Luogu https://www.luogu.com.cn/problem/P3338

---

Description

\[\forall i=1,2,\cdots,n \]

，求：

\[E_i~=~\sum\limits_{j = 1}^{i - 1} \frac{q_j}{(i - j)^2}~-~\sum\limits_{j = i + 1}^{n} \frac{q_j}{(i - j)^2} \]

其中

\[n~\in~\mathbb{N}~\cap~[1,10^5], \]

\[\forall j~\in~\mathbb{N}~\cap~[1,n]~,~\quad q_i~\in~\mathbb{R}~\cap~(0,10^9) \]

---

Solution

怎麼說呢，

\[\sum_{j = 1}^{i - 1} \frac{q_j}{(i - j)^2} \]

算是個卷積吧
那右邊那個反正直覺就是可以變成卷積
雖然但是，不難發現把 \(i\) 改成倒過來的馬上就得到一個卷積了
（比如說，原來的 \(i=n\) 就變成 \(i=1\)）

更加具體地，左邊是
首先設 \(f(j)=q_j\) 以及 \(g(x)=x^{-2}\)

\[\sum\limits_{j=1}^{i-1}f(j)g(i-j) \]

右邊是，

\[\sum\limits_{t=1}^{n-i}\frac{q_{i+t}}{t^2} \]

不妨設

\[h(x)\quad s.t.\quad h(n-i-t)=q_{i+t} \]

……至於下標不夠怎麼處理？
補0就可以啦（）

一度有過這樣的神奇想法
“為什麼輸出答案的時候，左邊和右邊兩個卷積不用先求和呢？”
……沒救了。卷積啊

---

經典錯誤了
就是那個
lim
。。因為卷積所以至少不小於 \(2n\)
然後因為
for(0; <lim; )
所以實際上 lim 不能比 \(2n+1\) 小。

---

Code

#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<cctype>
#include<cmath>

using namespace std;

const int MAXN = 4e5 + 10;
const double tpi = acos(-1.0);

int n, lim = 1, loglim = 0;
int rev[MAXN];

struct complex
{
    double real, imag;

    complex(const double &aug1 = 0, const double &aug2 = 0)
    {
        real = aug1;
        imag = aug2;
    }
    
    complex operator * (complex b)
    {
        return complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real);
    }

    complex operator + (complex b)
    {
        return complex(real + b.real, imag + b.imag);
    }
    
    complex operator - (complex b)
    {
        return complex(real - b.real, imag - b.imag);
    }

    void operator *= (complex b)
    {
        *this = *this * b;
    }

    void operator += (complex b)
    {
        real += b.real;
        imag += b.imag;
    }

    void operator -= (complex b)
    {
        real -= b.real;
        imag -= b.imag;
    }

    void operator *= (const double &x)
    {
        real *= x;
        imag *= x;
    }

} f[MAXN], g[MAXN], h[MAXN], Wn[MAXN][2];

void fft(complex *f, int opt)
{
    static complex t;

    for (int i = 0; i < lim; ++i)
    {
        if (rev[i] > i) swap(f[i], f[rev[i]]);
    }

    for (int n, stats, half_n = 1; half_n < lim; half_n = n)
    {
        n = half_n << 1;
        stats = lim / n;
        for (int pos = 0; pos < lim; pos += n)
        {
            for (int k = 0; k < half_n; ++k)
            {
                t = Wn[stats * k][opt] * f[pos + half_n + k];
                f[pos + half_n + k] = f[pos + k] - t;
                f[pos + k] += t;
            }
        }
    }

    if (opt == 1)
    {
        for (int i = 0; i < lim; ++i)
            f[i].real /= lim;
    }
}

void init()
{
    scanf("%d", &n);

    for (int i = 1; i <= n; ++i)
    {
        scanf("%lf", &f[i].real);
        h[n - i].real = f[i].real;
        g[i].real = 1.0 / i / i;
    }
    
    while (lim <= (n << 1))
    {
        lim <<= 1;
        ++loglim;
    }

    for (int i = 0; i < lim; ++i)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (loglim - 1));
    }

    const double modpi = 2.0 * tpi / lim;

    double theta;
    for (int i = 0; i < lim; ++i)
    {
        theta = modpi * i;
        Wn[i][0].real = cos(theta);
        Wn[i][1].real = +Wn[i][0].real;
        Wn[i][0].imag = sin(theta);
        Wn[i][1].imag = -Wn[i][0].imag;
    }
}

void work()
{
    fft(f, 0);
    fft(g, 0);
    fft(h, 0);

    for (int i = 0; i < lim; ++i)
    {
        f[i] *= g[i];
        h[i] *= g[i];
    }

    fft(f, 1);
    fft(h, 1);
}

int main()
{
    init();

    work();

    for (int i = 1; i <= n; ++i)
    {
        printf("%.3f\n", f[i].real - h[n - i].real);
    }

    system("pause");

    return 0;
}