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; }