P3338 [ZJOI2014]力（FFT）

P3338 [ZJOI2014]力

$$\begin{gathered} F_j=\sum _{i=1}^{j-1}\frac{q_i\times q_j}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i\times q_j}{(i-j{)}^2} \\ {E}_j=\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i}{(i-j{)}^2} \\ f(i)=q_i,g(i)=\frac{1}{i^2},f(0)=0,g(0)=0 \\ {E}_j=\sum _{i=0}^{j}f(i)g(j-i)-\sum _{i=j}^{n}f(i)g(i-j) \end{gathered}$$

前项是标准的多项式卷积，后项我们另考虑后项。

$\sum _{i=j}^{n}f(i)g(i-j)$，$T=i-j$，$i=T+j$，原式子$\sum _{i=0}^{n-j}f(i+j)g(i)$。

我们翻转$f(i)$得到$f^{\prime }(n-i)=f(i)$，代入原式子$\sum _{i=0}^{n-j}{f}^{\prime }(n-(i+j))g(i)$。

$n-j=T$，$\sum _{i=0}^T{f}^{\prime }(T-i)g(i)$，也就是一个多项式卷积了。

#include <bits/stdc++.h>

using namespace std;

struct Complex {
  double r, i;

  Complex(double _r = 0, double _i = 0) : r(_r), i(_i) {}
};

Complex operator + (const Complex &a, const Complex &b) {
  return Complex(a.r + b.r, a.i + b.i);
}

Complex operator - (const Complex &a, const Complex &b) {
  return Complex(a.r - b.r, a.i - b.i);
}

Complex operator * (const Complex &a, const Complex &b) {
  return Complex(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}

Complex operator / (const Complex &a, const Complex &b) {
  return Complex((a.r * b.r + a.i * b.i) / (b.r * b.r + b.i * b.i), (a.i * b.r - a.r * b.i) / (b.r * b.r + b.i * b.i));
}

typedef long long ll;

const int N = 1e6 + 10;

int r[N];

Complex a[N], b[N], c[N];

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void FFT(Complex *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  const double pi = acos(-1.0);
  for (int mid = 1; mid < lim; mid <<= 1) {
    Complex wn = Complex(cos(pi / mid), rev * sin(pi / mid));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      Complex w = Complex(1, 0);
      for (int k = 0; k < mid; k++, w = w * wn) {
        Complex x = f[cur + k], y = w * f[cur + mid + k];
        f[cur + k] = x + y, f[cur + mid + k] = x - y;
      }
    }
  }
  if (rev == -1) {
    for (int i = 0; i < lim; i++) {
      f[i].r /= lim;
    }
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  int n, lim = 1;
  scanf("%d", &n);
  for (int i = 1; i <= n; i++) {
    scanf("%lf", &a[i].r);
    b[n - i].r = a[i].r;
    c[i].r = 1.0 / i / i;
  }
  n <<= 1;
  while (lim <= n) {
    lim <<= 1;
  }
  get_r(lim);
  FFT(a, lim, 1);
  FFT(b, lim, 1);
  FFT(c, lim, 1);
  for (int i = 0; i < lim; i++) {
    a[i] = a[i] * c[i];
    b[i] = b[i] * c[i];
  }
  FFT(a, lim, -1);
  FFT(b, lim, -1);
  n >>= 1;
  for (int i = 1; i <= n; i++) {
    printf("%.3f\n", a[i].r - b[n - i].r);
  }
  return 0;
}