Problem
Description
给出 (n) 个数 (q_i),给出 (F_j) 的定义如下:
[F_j=sum_{i<j} frac{q_iq_j}{(i-j)^2} - sum_{i>j} frac{q_iq_j}{(i-j)^2}
]
令 (E_i=F_i/q_i),求 (E_i)。
Input Format
第一行一个整数(n)。
接下来 (n) 行每行输入一个数,第 (i) 行表示 (q_i)。
Output Format
(n) 行,第 (i) 行输出 (E_i)。与标准答案误差不超过 (10^{-2}) 即可。
Sample
Input
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
Output
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
Range
对于所有的数据,(nleq 100000, 0<q_i<10^9)。
Algorithm
多项式
Mentality
[E_i=sum_{j=1}^{i-1}frac{q_j}{(i-j)^2}-sum_{j=i+1}^nfrac{q_j}{(j-i)^2}
]
设 (g_i=frac{1}{i^2}) ,则有:
[E_i=sum_{j=1}^{i-1}q_jg_{i-j}-sum_{j=i+1}^nq_jg_{j-i}
]
单独算每个 (E_i) 的前一部分,再算所有 (E_i) 的后一部分。
由于 (q_0=g_0=0) ,则对于:
[sum_{j=1}^{i-1}q_jg_{i-j}=sum_{j=0}^{i}q_jg_{i-j}
]
直接 (FFT) 求得。
对于后一部分,将 (q) 数组翻转得到数组 (p) 。
则有:
[sum_{j=i+1}^nq_jg_{j-i}=sum_{j=0}^{n-i+1}p_jg_{n-i+1-j}
]
同样直接 (FFT) 计算。
Code
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <vector>
using namespace std;
#define cp complex<double>
const int Max_n = 3e5 + 5;
const double pi = acos(-1);
int n;
int rev[Max_n];
double dv[Max_n], ans[Max_n], q[Max_n];
cp f[Max_n], g[Max_n];
namespace FFT {
int lim, bit;
void dft(cp *f, int t) {
for (int i = 0; i < lim; i++)
if (rev[i] > i) swap(f[rev[i]], f[i]);
for (int len = 1; len < lim; len <<= 1) {
cp Wn = exp(cp(0, t * pi / len));
for (int i = 0; i < lim; i += len << 1) {
cp Wnk(1, 0);
for (int k = i; k < i + len; k++, Wnk *= Wn) {
cp x = f[k], y = Wnk * f[k + len];
f[k] = x + y, f[k + len] = x - y;
}
}
}
}
void fft(double *a, double *b, int tot) {
lim = 1, bit = 0;
while (lim <= tot) lim <<= 1, bit++;
for (int i = 0; i < lim; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
for (int i = 0; i < lim; i++) f[i] = a[i], g[i] = b[i];
dft(f, 1), dft(g, 1);
for (int i = 0; i < lim; i++) f[i] *= g[i];
dft(f, -1);
for (int i = 0; i < lim; i++) f[i] /= lim;
}
} // namespace FFT
using namespace FFT;
int main() {
#ifndef ONLINE_JUDGE
freopen("3338.in", "r", stdin);
freopen("3338.out", "w", stdout);
#endif
cin >> n;
for (int i = 1; i <= n; i++) scanf("%lf", &q[i]);
for (int i = 1; i <= n; i++) dv[i] = (double)(1.0 / i / i);
fft(q, dv, n << 1);
for (int i = 1; i <= n; i++) ans[i] += f[i].real();
reverse(q + 1, q + n + 1);
fft(q, dv, n << 1);
for (int i = 1; i <= n; i++)
printf("%.2lf
", (ans[i] -= f[n - i + 1].real()));
}