现在才明白卷积真的是the deep, dark fantasies……（逃

首先约去\(q_i\)，得到：

\[E_j = \sum_{i < j}\frac{q_i}{(j - i)^2} - \sum_{j < i}\frac{q_i}{(j - i)^2}\]

注意到如果很容易得到等式前半部分的高效求法，后半部分把序列反过来就能做了。

那么我们会发现，设\(g_x = \frac{1}{x^2}\)，然后式子前半部分（姑且称作\(p_j\)）可表示为：

\[p_j = \sum_{i < j} q_i g_{j - i}\]

这不就是个卷积吗？FFT一波即可。

代码：

/**************************************************************
    Problem: 3527
    User: danihao123
    Language: C++
    Result: Accepted
    Time:9984 ms
    Memory:28952 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <utility>
#include <queue>
#include <cassert>
#include <complex>
typedef long double R;
typedef std::complex<R> C;
const R eps = 1e-3;
int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x < 0) {
      return -1;
    } else {
      return 1;
    }
  }
}
const int maxn = 400005;
const double pi = acos(-1);
int bi, len;
int flip(int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if(x & (1 << i)) {
      ans += 1 << (bi - i - 1);
    }
  }
  return ans;
}
void fft(C *A, double g = 1) {
  for(int i = 0; i < len; i ++) {
    int R = flip(i);
#ifdef LOCAL
    // printf("The flipping of %d is %d\n", i, R);
#endif
    if(i < R) {
      std::swap(A[R], A[i]);
    }
  }
  for(int L = 1; L < len; L <<= 1) {
    C xi_n(cos(pi / (double(L))), sin(g * pi / (double(L))));
    for(int i = 0; i < len; i += L << 1) {
      C xi(1, 0);
      for(int j = i; j < i + L; j ++) {
        C a = A[j], b = A[j + L];
        A[j] = a + xi * b;
        A[j + L] = a - xi * b;
        xi = xi * xi_n;
      }
    }
  }
}
 
int main() {
  int n;
  static C A[maxn], rd[maxn]; static R ans[maxn];
  static R q[maxn];
  scanf("%d", &n);
  bi = 0; len = 1;
  while(len <= 2 * n) {
    len <<= 1;
    bi ++;
  }
  for(int i = 0; i < n; i ++) {
    scanf("%Lf", &q[i]); A[i] = q[i];
  }
  assert(sign(rd[0].real()) == 0);
  for(int i = 1; i < n; i ++) {
    rd[i] = 1.00 / ((R(i)) * (R(i)));
#ifdef LOCAL
    printf("rd[%d] : %.5Lf\n", i, rd[i].real());
#endif
  }
  fft(A); fft(rd);
  for(int i = 0; i < len; i ++) A[i] *= rd[i];
  fft(A, -1);
  for(int i = 0; i < n; i ++) {
    ans[i] += A[i].real() / len;
#ifdef LOCAL
    printf("delta_v of %d : %.5Lf\n", i, A[i].real() / len);
#endif
  }
  std::reverse(q, q + n);
  for(int i = 0; i < len; i ++) A[i] = 0;
  for(int i = 0; i < n; i ++) {
    A[i] = q[i];
  }
  fft(A);
  for(int i = 0; i < len; i ++) A[i] *= rd[i];
  fft(A, -1);
  for(int i = 0; i < n; i ++) {
    ans[i] -= A[n - 1 - i].real() / len;
    printf("%.3Lf\n", ans[i]);
  }
  return 0;
}

• 相关文章
• 全局相关文章