BZOJ3527 [Zjoi2014]力
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=\frac{F_i}{q_i}$,求$E_i$.
Input
第一行一个整数n。
接下来n行每行输入一个数,第i行表示qi。
n≤100000,0<qi<1000000000
Output
n行,第i行输出Ei。与标准答案误差不超过1e-2即可。
Sample Input
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
Sample Output
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
3439.793
7509018.566
4595686.886
10903040.872
题解
由于
$$E_i = \sum_{j<i}\frac{q_j}{(i-j)^2} - \sum_{j>i}\frac{q_j}{(i-j)^2}$$
我们令$t_i = [i>0]i^{-2}$,那么
$$E_i = \sum_{j<i}q_jt_{i-j} - \sum_{j>i}q_jt_{j-i}$$
对于前一项,直接FFT计算;对于后一项,可以将数组q翻转后即可得到卷积形式,FFT即可。
代码:
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdio>
typedef std::complex<double> Complex;
const double pi = acos(-1.0);
const int N = 400000;
double E[N], Q[N];
Complex A[N], B[N];
void FFT(Complex *P, int len, int opt) {
for (int i = 1, j = len >> 1; i < len; ++i) {
if (i < j) std::swap(P[i], P[j]);
int k = len;
do {
k >>= 1;
j ^= k;
} while (~j & k);
}
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * pi * opt / h), sin(2 * pi * opt / h));
for (int j = 0; j < len; j += h) {
Complex w(1.0, 0.0);
for (int k = 0; k < h / 2; ++k) {
Complex t1 = P[j + k] , t2 = P[j + k + h / 2] * w;
P[j + k] = t1 + t2;
P[j + k + h / 2] = t1 - t2;
w *= wn;
}
}
}
if (opt == -1) {
for (int i = 0; i < len; ++i)
P[i] /= len;
}
}
int main() {
int n;
scanf("%d", &n);
for (int i = 0; i < n; ++i) scanf("%lf", &Q[i]);
int len = 1;
while (len < 2 * n) len <<= 1;
for (int i = 0; i < n; ++i) A[i] = Q[i];
for (int i = n; i < len; ++i) A[i] = 0.0;
B[0] = 0.0;
for (int i = 1; i < n; ++i) B[i] = 1.0 / i / i;
for (int i = n; i < len; ++i) B[i] = 0.0;
FFT(A, len, 1);
FFT(B, len, 1);
for (int i = 0; i < len; ++i) A[i] *= B[i];
FFT(A, len, -1);
for (int i = 0; i < n; ++i) E[i] += A[i].real();
for (int i = 0; i < n; ++i) A[i] = Q[n - i - 1];
for (int i = n; i < len; ++i) A[i] = 0.0;
B[0] = 0.0;
for (int i = 1; i < n; ++i) B[i] = 1.0 / i / i;
for (int i = n; i < len; ++i) B[i] = 0.0;
FFT(A, len, 1);
FFT(B, len, 1);
for (int i = 0; i < len; ++i) A[i] *= B[i];
FFT(A, len, -1);
for (int i = 0; i < n; ++i) E[n - i - 1] -= A[i].real();
for (int i = 0; i < n; ++i) printf("%.10lf\n", E[i]);
return 0;
}

浙公网安备 33010602011771号