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

Sample Output

-16838672.693
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
", E[i]);
  return 0;
}