洛谷 [P3338] 力

FFT

\[E_i = F_i / q_i = \sum_{i<j} \frac {q_j} {(i - j)^2} - \sum _{ i > j} \frac{q _ j} {(i - j)^2} \]

设 \(p _ i = q_{n - i}\) \(g(i) = \frac 1 {i^2}\)

则 \(E_i = \sum_{j=1}^{i-1} q _ i g(j - i) - \sum _ {j = i + 1} ^ n q_j g(i - j)\)

即 \(E_i = \sum_{j=1}^{i-1} q_i g(i-1) - \sum_{j=1}^{n - i} p_i g(i - j)\)

我们发现这是一个卷积的形式, 可以用 FFT 处理

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int MAXN = 400005;
const double PI = acos(-1);
struct Complex {
	double x, y;
	Complex(double xx = 0.0, double yy = 0.0) {
		x = xx; y = yy;
	}
	Complex operator + (const Complex &u) const {
		return Complex(x + u.x, y + u.y);
	}
	Complex operator - (const Complex &u) const {
		return Complex(x - u.x, y - u.y);
	}
	Complex operator * (const Complex &u) const{
		return Complex(x * u.x - y * u.y, x * u.y + y * u.x);
	}
}a[MAXN], b[MAXN], c[MAXN], buf[MAXN];
int n, lim = 1, limcnt, rev[MAXN];
void fft(Complex a[], int opt) {
	for(int i = 0; i <= lim; i++) {
		if(i < rev[i]) swap(a[i], a[rev[i]]);
	}
	for(int mid = 1; mid < lim; mid <<= 1) {
		Complex wn = Complex(cos(PI / mid), opt * sin(PI / mid));
		for(int R = mid << 1, j = 0; j < lim; j += R) {
			Complex w = Complex(1.0, 0.0);
			for(int k = 0; k < mid; k++) {
				Complex x = a[j + k], y = w * a[j + mid + k];
				a[j + k] = x + y;
				a[j + mid + k] = x - y;
				w = w * wn;
			}
		}
	}
	if(opt == -1) {
		for(int i = 0; i <= lim; i++) a[i].x = a[i].x / lim;
	}
}
int main() {
	cin >> n;
	for(int i = 1; i <= n; i++) {
		scanf("%lf", &a[i].x);
		b[n - i].x = a[i].x;
		c[i].x = (1.0 / (double)i) / (double)i;
	}
	while(lim <= (n << 1)) lim <<= 1, limcnt++;
 	for(int i = 0; i <= lim; i++) {
 		rev[i] = (rev[i>>1]>>1) | ((i&1)<<(limcnt - 1));
 	}
 	fft(a, 1);fft(b, 1); fft(c, 1);
 	for(int i = 0; i <= lim; i++) a[i] = a[i] * c[i];
 	for(int i = 0; i <= lim; i++) b[i] = b[i] * c[i];
 	fft(a, -1); fft(b, -1);
 	for(int i = 1; i <= n; i++) {
 		printf("%.4f\n", a[i].x - b[n - i].x);
 	}
	return 0;
}```