Description

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

，求：

[E_i~=~sumlimits_{j = 1}^{i - 1} frac{q_j}{(i - j)^2}~-~sumlimits_{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})

[sumlimits_{j=1}^{i-1}f(j)g(i-j) ]

右边是，

[sumlimits_{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
", f[i].real - h[n - i].real);
    }

    system("pause");

    return 0;
}