Luogu3338 [ZJOI2014]力

题目蓝链

Description

给定一个长度为(n)的数列(q_i)，定义

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

求出所有的(E_i)

Solution

直接上推导过程

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

于是我们可以令

[F(x) = q_x \ G(x) = egin{cases} -frac{1}{x^2}~~~(x < 0) \ 0~~~~~~~~(x = 0) \ frac{1}{x^2}~~~~~~(x > 0) end{cases} ]

则(下标默认从0开始)

[E_j = sum_{i = 0}^{n - 1} F(j) cdot G(j - i) ]

所以这显然是一个卷积形式，所以我们就直接FFT求出卷积的系数即(E_i)即可

Code

#include <bits/stdc++.h>

using namespace std;

#define fst first
#define snd second
#define mp make_pair
#define squ(x) ((LL)(x) * (x))
#define debug(...) fprintf(stderr, __VA_ARGS__)

typedef long long LL;
typedef pair<int, int> pii;

template<typename T> inline bool chkmax(T &a, const T &b) { return a < b ? a = b, 1 : 0; }
template<typename T> inline bool chkmin(T &a, const T &b) { return a > b ? a = b, 1 : 0; }

inline int read() {
	int sum = 0, fg = 1; char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') fg = -1;
	for (; isdigit(c); c = getchar()) sum = (sum << 3) + (sum << 1) + (c ^ 0x30);
	return fg * sum;
}

namespace FFT {

	const int MAX_LEN = 1 << 19;
	const double PI = acos(-1.0);

	struct com {
		double a, b;
		com (double _a = 0.0, double _b = 0.0): a(_a), b(_b) { }
		com operator + (const com &t) const { return com(a + t.a, b + t.b); }
		com operator - (const com &t) const { return com(a - t.a, b - t.b); }
		com operator * (const com &t) const { return com(a * t.a - b * t.b, a * t.b + b * t.a); }
	};

	int cnt, len, rev[MAX_LEN];
	com g[MAX_LEN];

	void init(int N) {
		for (cnt = -1, len = 1; len <= N; len <<= 1) ++cnt;
		for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << cnt);
		g[0] = com(1.0, 0.0);
		com G(cos(PI * 2 / len), sin(PI * 2 / len));
		for (int i = 1; i < len; i++) g[i] = g[i - 1] * G;
	}

	void DFT(com *x, int op) {
		for (int i = 0; i < len; i++) if (i < rev[i]) swap(x[i], x[rev[i]]);
		for (int k = 2; k <= len; k <<= 1)
			for (int j = 0; j < len; j += k)
				for (int i = 0; i < k / 2; i++) {
					com X = x[j + i], Y = x[j + i + k / 2] * g[len / k * (~op ? i : (i ? k - i : i))];
					x[j + i] = X + Y, x[j + i + k / 2] = X - Y;
				}
		if (op < 0) for (int i = 0; i < len; i++) x[i].a /= len;
	}

	void mul(double *a, int n, double *b, int m, double *c) {
		init(n + m);
		static com F[MAX_LEN], G[MAX_LEN], S[MAX_LEN];
		for (int i = 0; i < len; i++) F[i] = com(i <= n ? a[i] : 0.0, 0.0);
		for (int i = 0; i < len; i++) G[i] = com(i <= m ? b[i] : 0.0, 0.0);
		DFT(F, 1), DFT(G, 1);
		for (int i = 0; i < len; i++) S[i] = F[i] * G[i];
		DFT(S, -1);
		for (int i = 0; i <= n + m; i++) c[i] = S[i].a;
	}

}

const int maxn = 3e5 + 10;

int n;
double f[maxn], g[maxn], e[maxn];

int main() {
#ifdef xunzhen
	freopen("force.in", "r", stdin);
	freopen("force.out", "w", stdout);
#endif

	n = read();
	for (int i = 0; i < n; i++) scanf("%lf", &f[i]);
	for (int i = 0; i < (n << 1); i++) if (i != n) g[i] = (i < n ? -1.0 : 1.0) / squ(i - n);

	FFT::mul(f, n, g, (n << 1) - 1, e);

	for (int i = n; i < (n << 1); i++) printf("%.3lf
", e[i]);

	return 0;
}