[ZJOI2014]力 FFT

题面

题解：

\[F_j = \sum_{i < j}\frac{q_iq_j}{(i – j)^2} – \sum_{i > j}{\frac{q_iq_j}{(i – j)^2}}
\]

\[E_j = \sum_{i < j}\frac{q_i}{(i – j)^2} – \sum_{i > j}{\frac{q_i}{(i – j)^2}}
\]

对式子的2个部分分别计算。

令\(S_i = i^2\)

\[\sum_{i < j}\frac{q_i}{(i – j)^2} = \sum_{i < j}q_i S_{j – i}
\]

看上去就是卷积形式，FFT计算即可。

对于后半部分，将序列翻转，\(i > j\)就变成\(i < j\)了，而\(S\)可以看做距离，所以不会变，直接计算就好了.

计算完之后需要将序列翻转回来

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define ld double
#define LL long long
#define AC 310000const double pi = acos(-1);
int n, lim = 1, len;
int Next[AC];
ld q[AC], f[AC];struct node{
    ld x, y;
    node(ld xx = 0, ld yy = 0){x = xx, y = yy;}
}a[AC], b[AC], s[AC];node operator * (node x, node y){return node(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);}
node operator + (node x, node y){return node(x.x + y.x, x.y + y.y);}
node operator - (node x, node y){return node(x.x - y.x, x.y - y.y);}inline int read()
{
    int x = 0;char c = getchar();
    while(c > '9' || c < '0') c = getchar();
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x;
}void FFT(node *A, int opt)
{
    for(R i = 0; i < lim; i ++)
        if(i < Next[i]) swap(A[i], A[Next[i]]);
    for(R i = 1; i < lim ; i <<= 1)
    {
        node W(cos(pi / i), opt * sin(pi / i));
        for(R r = i << 1, j = 0; j < lim; j += r)
        {
            node w(1, 0);
            for(R k = 0; k < i ; k ++, w = w * W)
            {
                node x = A[j + k], y = w * A[j + k + i];
                A[j + k] = x + y, A[j + k + i] = x - y;
            }
        }
    }
}void pre()
{
    n = read() - 1;
    for(R i = 0; i <= n; i ++) scanf("%lf", &q[i]);
    while(lim <= n + n) lim <<= 1, ++ len;
    for(R i = 0; i <= lim; i ++)
        Next[i] = (Next[i >> 1] >> 1) | ((i & 1) << (len - 1));
}void work()
{
    for(R i = 0; i <= n; i ++)
    {
        if(i != 0) s[i].x = 1.0 / i / i;
        a[i].x = q[i], b[n - i].x = q[i];
    }
    FFT(s, 1);
    FFT(a, 1);
    for(R i = 0; i < lim; i ++) a[i] = a[i] * s[i];
    FFT(a, -1);
    for(R i = 0; i < lim; i ++) f[i] = a[i].x / lim;
    FFT(b, 1);
    for(R i = 0; i < lim; i ++) b[i] = b[i] * s[i];
    FFT(b, -1);
    for(R i = 0; i <= n; i ++)
        if(i < n - i) swap(b[i], b[n - i]);
    for(R i = 0; i < lim; i ++) f[i] -= b[i].x / lim;
    for(R i = 0; i <= n; i ++) printf("%.3lf\n", f[i]);
}int main()
{
    //freopen("in.in", "r", stdin);
    pre();
    work();
    //fclose(stdin);
    return 0;
}