首先我们知道(displaystyle E_j=sum_{i<j}frac{q_i}{(i-j)^2}-sum_{i>j}frac{q_i}{(i-j)^2}),
设(displaystyle g[i]=frac{1}{i^2}),因为(g)是偶函数，所以(displaystyle E_j=sum_{i=0}^{j-1} q_i g[j-i]-sum_{i=j+1}^n q_ig[j-i])。
前面这东西很明显就是卷积，处理后面就要发挥人类智慧了，将(q_i) 翻转成新数组(q_i')，原来的式子就成了
(displaystyle E_j=sum_{i=0}^{j-1} q_i g[j-i]-sum_{i=0}^{j-1} q_i'g[j-i]),
后面这东西也变成卷积啦，用FFT处理一下就ok啦。
时间复杂度O(n log n).
#include<iostream>
#include<cstdio>
#include<cmath>
#define DB double
using namespace std;
int n, lim = 1, tmp;
const DB PI = acos(-1);
const int N = 400010;
struct lmaginary
{
DB x, y;
lmaginary(double X = 0, double Y = 0) {x = X, y = Y;}
friend lmaginary operator +(const lmaginary &a, const lmaginary &b)
{return (lmaginary) {a.x + b.x, a.y + b.y};}
friend lmaginary operator -(const lmaginary &a, const lmaginary &b)
{return (lmaginary) {a.x - b.x, a.y - b.y};}
friend lmaginary operator *(const lmaginary &a, const lmaginary &b)
{return (lmaginary) {a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x};}
} f1[N], f2[N], g[N];
DB q[N];
int r[N];
void FFT(lmaginary *A, int lim, int opt)
{
for (int i = 0; i < lim; ++i)
if (i < r[i])swap(A[i], A[r[i]]);
for (int mid = 1; mid < lim; mid <<= 1)
{ //长度的一半
lmaginary wn(cos(PI / mid), opt * sin(PI / mid));
for (int len = mid << 1, j = 0; j < lim; j += len)
{
lmaginary w((DB)1, (DB)0);
for (int k = j; k < mid + j; ++k, w = w * wn)
{
lmaginary a = A[k], b = w * A[k + mid];
A[k] = a + b;
A[k + mid] = a - b;
}
}
}
}
void YYCH()
{
for (int i = 1; i <= n; ++i)
{
g[i].x = (1.0 / i / i);
f1[i].x = f2[n - i + 1].x = q[i];
}
while (lim < 2 * n)lim <<= 1, ++tmp;
for (int i = 0; i < lim; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (tmp - 1));
}
int main()
{
cin >> n;
for (int i = 1; i <= n; ++i)scanf("%lf", &q[i]);
YYCH();
FFT(f1, lim, 1); FFT(f2, lim, 1); FFT(g, lim, 1);
for (int i = 0; i < lim; ++i)
f1[i] = f1[i] * g[i], f2[i] = f2[i] * g[i];
FFT(f1, lim, -1); FFT(f2, lim, -1);
for (int i = 1; i <= n; ++i)
printf("%f
", (f1[i].x - f2[n - i + 1].x) / lim);
return 0;
}