题目链接
题解
首先我们需要对这个式子进行化简,否则对着这么大一坨东西只能暴力。。。
\[F_i=\sum_{j<i} \frac{q_iq_j}{(i-j)^2}-\sum_{j>i} \frac{q_iq_j}{(i-j)^2}\]
根据题目给出的定义,带入\(E\)中
\[E_i=\sum_{j=1}^{i-1}\frac{q_j}{(i-j)^2}-\sum_{j=i+1}^{n}\frac{q_j}{(j-i)^2}\]
形式稍微改变了一下,本质一样
需要处理这个分母上的两个不好看的小东西
设\(f_i=\frac{1}{i^2}\)
那么可以得到原式变成
\[E_i=\sum^{i-1}_{j=1}q_jf_{i-j}-\sum^n_{j=i+1}q_jf_{j-i}\]
将这个式子拆成两部分
先讨论\(-\)号前面一部分
可以发现,减号前一部分就是一个卷积的形式,那么直接FFT
套上去就可以了。【FFT学习笔记(万字大篇)】
对于后一部分,认真观察可以发现
形式和前一部分非常相似,其实就是相反了一下。
我们将原来的数组倒序存储一下,和前一部分做一样的FFT
就可以了。
代码
#include <bits/stdc++.h>
#define ms(a, b) memset(a, b, sizeof(a))
#define ll long long
#define ull unsigned long long
#define ms(a, b) memset(a, b, sizeof(a))
#define inf 0x3f3f3f3f
#define db double
#define Pi acos(-1)
#define eps 1e-8
#define N 400005
using namespace std;
template <typename T> T sqr(T x) { return x * x; }
template <typename T> void read(T &x) {
x = 0; T fl = 1; char ch = 0;
for (; ch < '0' || ch > '9'; ch = getchar()) if (ch == '-') fl = -1;
for (; ch >= '0' && ch <= '9'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48);
x *= fl;
}
template <typename T> T power(T x, T y, T mod) { T res = 1; for (; y; y >>= 1) { if (y & 1) res = (res * x) % mod; x = (x * x) % mod; } return res; }
template <typename T> void write(T x) { if (x < 0) x = -x, putchar('-'); if (x > 9) write(x / 10); putchar(x % 10 + '0'); }
template <typename T> void writeln(T x) { write(x); puts(""); }
struct Complex {
db x, y;
Complex (db xx = 0, db yy = 0) { x = xx; y = yy; }
Complex operator + (Complex B) { return Complex(x + B.x, y + B.y); }
Complex operator - (Complex B) { return Complex(x - B.x, y - B.y); }
Complex operator * (Complex B) { return Complex(x * B.x - y * B.y, x * B.y + y * B.x); }
}a[N], b[N], c[N], d[N];
int r[N];
int n, m, limit, l;
void FFT(Complex *a, int type) {
for (int i = 0; i < limit; i ++) if (i < r[i]) swap(a[i], a[r[i]]);
for (int mid = 1; mid < limit; mid <<= 1) {
Complex wn(cos(Pi / mid), type * sin(Pi / mid));
for (int R = mid << 1, j = 0; j < limit; j += R) {
Complex w(1, 0);
for (int k = 0; k < mid; k ++, w = w * wn) {
Complex x = a[j + k], y = w * a[j + mid + k];
a[j + k] = x + y; a[j + mid + k] = x - y;
}
}
}
if (type == -1) for (int i = 0; i < limit; i ++) a[i].x /= 1.0 * limit;
}
int main() {
read(n);
for (int i = 1; i <= n; i ++) {
db x; scanf("%lf", &x);
a[i].x = c[n - i + 1].x = x; b[i].x = d[i].x = 1.0 / sqr(i * 1.0);
}
limit = 1; while (limit <= (n << 1)) limit <<= 1, l ++;
for (int i = 0; i < limit; i ++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
FFT(a, 1); FFT(b, 1);
for (int i = 0; i < limit; i ++) a[i] = a[i] * b[i];
FFT(a, -1);
FFT(c, 1); FFT(d, 1);
for (int i = 0; i < limit; i ++) c[i] = c[i] * d[i];
FFT(c, -1);
for (int i = 1; i <= n; i ++) printf("%.3lf\n", a[i].x - c[n - i + 1].x);
return 0;
}