题目描述
给出 $n$ 个数 $q_i$ ，给出 $F_j$ 的定义如下：
$$F_j = \sum_{i<j} \frac{q_i q_j}{{(i-j)}^2} – \sum_{i>j} \frac{q_i q_j}{{(i-j)}^2} $$
令 $E_i = \frac{F_i}{q_i} $，求 $E_i$
输入格式
第一行一个整数 $n$
接下来 $n$ 行每行输入一个整数，第 $i$ 行表示 $q_i$
输出格式
$n$ 行，第 $i$ 行输出 $E_i$
样例输入
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
样例输出
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
数据范围及说明
时间限制 : $3s$
空间限制 : $128MB$
对于 $30 \%$ 的数据，$n \leq 1000$
对于 $50 \%$ 的数据，$n \leq 60000$
对于 $100 \%$ 的数据，$n \leq 100000, 0 < q_i < 1000000000 $
有 $spj$，误差不超过 $0.01$ 即可
解题思路
$$
\begin{aligned}
E_i &= \frac{F_i}{q_i} \\
&= \frac{\sum_{j < i} \frac{q_iq_j}{{(i-j)}^2} – \sum_{j>i} \frac{q_iq_j}{{(i-j)}^2}}{q_i} \\
&= \sum_{j < i} \frac{q_j}{{(i-j)}^2} – \sum_{j>i} \frac{q_j}{{(i-j)}^2} \\
&= \sum_{j=1}^{i-1} \frac{q_j}{{(i-j)}^2} – \sum_{j=i+1}^{n} \frac{q_j}{{(i-j)}^2} \\
\end{aligned}
$$
令 $f[i] = q_i, g[i] = \frac{1}{i^2}$，令 $f[0] = 0, g[0] = 0$
$$
\begin{aligned}
E_i &= \sum_{j=1}^{i-1} f[j]g[i-j] – \sum_{j=i+1}^{n} f[j]g[j-i] \\
&= \sum_{j=1}^{i-1} f[j]g[i-j] – \sum_{j=1}^{n-i} f[j+i]g[j] \\
&= \sum_{j=0}^{i-1} f[j]g[i-j] – \sum_{j=0}^{n-i} f[j+i]g[j] \\
\end{aligned}
$$
翻转 $f$ 数组得到 $f’$
$$
\begin{aligned}
E_i &= \sum_{j=0}^{i-1} f[j] g[i-j] – \sum_{j=0}^{n-i} f'[n-i-j+1]g[j] \\
\end{aligned}
$$
式子转化为两个卷积的形式，$FFT$ 即可。
代码
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int MAXN = 262144 + 10;
const long double PI = acos(-1.0);
struct complex{
long double x, y;
complex(long double _x = 0, long double _y = 0){
x = _x, y = _y;
}
};
complex operator + (complex a, complex b){
return complex(a.x + b.x, a.y + b.y);
}
complex operator - (complex a, complex b){
return complex(a.x - b.x, a.y - b.y);
}
complex operator * (complex a, complex b){
return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}
complex omega[MAXN], omegaInverse[MAXN];
int place[MAXN];
inline void init(int n){
int k = 0;
while((1 << k) < n)
k++;
for(register int i=0; i<n; i++){
place[i] = (place[i >> 1] >> 1) | ((i & 1) << (k - 1));
}
long double single = 2 * PI / n;
for(register int i=0; i<n; i++){
omega[i] = complex(cos(single * i), sin(single * i));
omegaInverse[i] = complex(omega[i].x, -omega[i].y);
}
}
inline void transform(complex *a, int n, complex *omega){
for(register int i=0; i<n; i++){
if(i < place[i]){
swap(a[i], a[place[i]]);
}
}
for(register int range=2; range <= n; range <<= 1){
int mid = range >> 1;
int single = n / range;
for(register complex *data=a; data != a + n; data += range){
for(register int i=0; i<mid; i++){
complex tmp = omega[single * i] * data[i + mid];
data[i + mid] = data[i] - tmp;
data[i] = data[i] + tmp;
}
}
}
}
int n, m = 1;
long double f[MAXN], fInverse[MAXN], g[MAXN];
long double resultA[MAXN], resultB[MAXN];
complex tmp[MAXN], tmp2[MAXN];
inline void work(long double *a, long double *b, long double *result){
for(register int i=0; i<m; i++){
tmp[i].x = a[i];
tmp2[i].x = b[i];
tmp[i].y = 0;
tmp2[i].y = 0;
}
transform(tmp, m, omega);
transform(tmp2, m, omega);
for(register int i=0; i<m; i++){
tmp[i] = tmp[i] * tmp2[i];
}
transform(tmp, m, omegaInverse);
for(register int i=0; i<n; i++){
result[i] = tmp[i].x / m;
}
}
int main(){
scanf("%d", &n);
for(register int i=1; i<=n; i++){
scanf("%Lf", &f[i]);
g[i] = (long double)1 / i / i;
fInverse[i] = f[i];
}
n++;
while(m < (n << 1)){
m <<= 1;
}
init(m);
reverse(fInverse + 1, fInverse + n);
work(f, g, resultA);
work(fInverse, g, resultB);
n--;
for(register int i=1; i<=n; i++){
printf("%.10Lf\n", resultA[i] - resultB[n - i + 1]);
}
return 0;
}