两道FFT题目略解

[ZJOI2014]力

Description

给出 (n) 个数 (q_1,q_2, dots q_n)，定义

[F_j = sum_{i = 1}^{j - 1} frac{q_i imes q_j}{(i - j)^2} - sum_{i = j + 1}^{n} frac{q_i imes q_j}{(i - j)^2} ]

[E_i = frac{F_i}{q_i} ]

对 (1 leq i leq n)，求 (E_i) 的值。

Solution

[E_j = frac{F_j}{q_j} \ =sum_{i = 1}^{j - 1} frac{q_i}{(j - i)^2}~-~sum_{i = j + 1}^{n} frac{q_i}{(i - j)^2} ]

然后，我们令

[f_i = frac{1}{i!},g_i=q_i ]

设(g'_i)表示将(g_i)反向后的数组，那么原式就变成了

[E_j =sum_{i = 1}^{j - 1} f_{j-i}g_i - sum_{i = 1}^{n-j}f_{j-i}g'_i ]

很显然这是一个卷积的形式，直接(FFT)求即可。

Code

#include <bits/stdc++.h>
using namespace std;

inline int ty() {
	char ch = getchar(); int x = 0, f = 1;
	while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
	while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
	return x * f;
}

const int _ = 4e5 + 10;
const double Pi = acos(-1.0);
struct Complex {
	double x, y;
	Complex operator+(const Complex &b) const { return {x + b.x, y + b.y}; }
	Complex operator-(const Complex &b) const { return {x - b.x, y - b.y}; }
	Complex operator*(const Complex &b) const { return {x * b.x - y * b.y, x * b.y + y * b.x}; }
} F[_], G[_];
int N, lim = 1, pos[_];
double p[_], q[_], t[_], a[_], b[_];

void fft(int lim, Complex *a, int op) {
	for (int i = 0; i < lim; ++i)
		if (i < pos[i]) swap(a[i], a[pos[i]]);
	for (int len = 2; len <= lim; len <<= 1) {
		int mid = len >> 1;
		Complex Wn = { cos(2.0 * Pi / len), op * sin(2.0 * Pi / len) };
		for (int i = 0; i < lim; i += len) {
			Complex w = { 1, 0 };
			for (int j = 0; j < mid; ++j, w = w * Wn) {
				Complex x = a[i + j], y = w * a[i + j + mid];
				a[i + j] = x + y;
				a[i + j + mid] = x - y;
			}
		}
	}
}

void work(double *a, double *b, double *ret) {
	for (int i = 0; i < lim; ++i) {
		F[i].x = a[i], G[i].x = b[i];
		F[i].y = G[i].y = 0;
	}
	fft(lim, F, 1);
	fft(lim, G, 1);
	for (int i = 0; i < lim; ++i) F[i] = F[i] * G[i];
	fft(lim, F, -1);
	for (int i = 1; i <= N; ++i) ret[i] = F[i].x / lim;
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("force.in", "r", stdin);
	freopen("force.out", "w", stdout);
#endif
	scanf("%d", &N);
	for (int i = 1; i <= N; ++i) {
		scanf("%lf", &p[i]);
		q[i] = p[i], t[i] = 1.0 / i / i;
	}
	reverse(q + 1, q + N + 1);
	int k = 0; while (lim <= N + N) lim <<= 1, ++k;
	for (int i = 0; i < lim; ++i) pos[i] = ((pos[i >> 1] >> 1)) | ((i & 1) << (k - 1));
	work(p, t, a);
	work(q, t, b);
	for (int i = 1; i <= N; ++i) printf("%.3lf
", a[i] - b[N - i + 1]);
	return 0;
}

[AH2017/HNOI2017]礼物

Description

传送门

Solution

考虑给数列加上一个增量后，将式子进行展开：

[sum_{i=1}^{n} (a_i-b_i+x)^2 = \ sum_{i=1}^{n}a_i^2 + sum_{i=1}^{n}b_i^2 + nx^2 + 2x (sum_{i=1}^n a_i - sum_{i=1}^{n}b_i) - 2sum_{i=1}^{n} a_i b_i ]

前面的项最小值是确定的，直接利用二次函数求解即可。

现在问题就变成了求(sum_{i=1}^{n} a_i b_i)的最大值。

另外还有关于环的问题，所以只需将(b_i)反向，然后倍长(a_i)，求出卷积后，取第(n+1)项到第(2n)项的最大值即可。

Code

#include <bits/stdc++.h>
using namespace std;

namespace IO {
const int bufl = 1 << 15;
char buf[bufl], *s = buf, *t = buf;
inline int fetch() {
  if (s == t) {
    t = (s = buf) + fread(buf, 1, bufl, stdin);
    if (s == t) return EOF;
  }
  return *s++;
}
inline int ty() {
  int a = 0;
  int b = 1, c = fetch();
  while (!isdigit(c)) b ^= c == '-', c = fetch();
  while (isdigit(c)) a = a * 10 + c - 48, c = fetch();
  return b ? a : -a;
}
}  // namespace IO
using IO::ty;

typedef long long ll;
const int _ = 1e6 + 10;
const double Pi = acos(-1.0);
struct Complex {
	double x, y;
	Complex operator+(const Complex &rhs) const { return { x + rhs.x, y + rhs.y }; }
	Complex operator-(const Complex &rhs) const { return { x - rhs.x, y - rhs.y }; }
	Complex operator*(const Complex &rhs) const { return { x * rhs.x - y * rhs.y, x * rhs.y + y * rhs.x }; }
} a[_], b[_];
int N, M, MX, maxx, r[_];
ll ans, suma, sumb;

inline ll _2(ll x) { return x * x; }

void fft(int lim, Complex *a, int op) {
	for (int i = 0; i < lim; ++i)
		if (i < r[i]) swap(a[i], a[r[i]]);
	for (int len = 2; len <= lim; len <<= 1) {
		int mid = len >> 1;
		Complex Wn = { cos(Pi / mid), op * sin(Pi / mid) };
		for (int i = 0; i < lim; i += len) {
			Complex w = { 1, 0 };
			for (int j = 0; j < mid; ++j, w = w * Wn) {
				Complex x = a[i + j], y = w * a[i + j + mid];
				a[i + j] = x + y;
				a[i + j + mid] = x - y;
			}
		}
	}
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("gift.in", "r", stdin);
	freopen("gift.out", "w", stdout);
#endif
	N = ty(), MX = ty();
	for (int i = 1; i <= N; ++i) {
		a[i + N].x = a[i].x = ty();
		ans += _2((ll)a[i].x);
		suma += (ll)a[i].x;
	}
	for (int i = N; i >= 1; --i) {
		b[i].x = ty();
		ans += _2((ll)b[i].x);
		sumb += (ll)b[i].x;
	}
	M = N, N <<= 1;

	ll m1 = floor(1.0 * (sumb - suma) / N);
	ll m2 = ceil(1.0 * (sumb - suma) / N);
	ll t = suma - sumb;
	ans += min(2ll * m1 * t + M * m1 * m1, 2ll * m2 * t + M * m2 * m2);

	int lim = 1, k = 0;
	while (lim <= N + M) lim <<= 1, ++k;
	for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (k - 1));
	fft(lim, a, 1);
	fft(lim, b, 1);
	for (int i = 0; i < lim; ++i) a[i] = a[i] *  b[i];
	fft(lim, a, -1);
	for (int i = 0; i < lim; ++i) a[i].x = (ll)(a[i].x / lim + 0.5);

	ll maxx = 0;
	for (int i = M + 1; i <= M + M; ++i) maxx = max(maxx, (ll)a[i].x);
	ans -= 2 * maxx;
	printf("%lld
", ans);
	return 0;
}