题目：

BZOJ3527

分析：

FFT应用第一题……

首先很明显能把\(F_j\)约掉，变成：

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

然后去膜拜题解，我们知道两个多项式相乘的方式如下：

\[C_j=\sum_{i=0}^j A_iB_{j-i}
\]

那么，如果把\(E_j\)的表达式化成上面那个形式，就可以用FFT计算了。（不会FFT？戳我：【知识总结】快速傅里叶变换（FFT））

先看减号前面的部分。显然可以变成（为了叙述方便，读入的\(q\)的下标为\([0,n)\)）：

\[C_j=\sum_{i=0}^j F_iG_{j-i}
\]

其中\(F_i=q_i\)，\(G_i=\frac{1}{i^2}\)。特别地，\(G_0=0\)。

减号后面要处理\(j\)位置以后的，怎么办？大力把\(q\)数组翻过来，这样就相当于求\(n-j-1\)以前的了：

\[D_j=\sum_{i=0}^{j} F'_iG_{j-i}
\]

其中\(F'_j=q_{n-j-1}\)

那么答案\(E_j=C_j-D_{n-j-1}\)

代码：

注意\(g\)要初始化……

#include <cstdio>

#include <algorithm>

#include <cstring>

#include <cctype>

#include <cmath>

using namespace std;

#define _ 0

namespace zyt

{

	template<typename T>

	inline void read(T &x)

	{

		char c;

		bool f = false;

		x = 0;

		do

			c = getchar();

		while (c != '-' && !isdigit(c));

		if (c == '-')

			f = true, c = getchar();

		do

			x = x * 10 + c - '0', c = getchar();

		while (isdigit(c));

		if (f)

			x = -x;

	}

	inline void read(double &x)

	{

		scanf("%lf", &x);

	}

	template<typename T>

	inline void write(T x)

	{

		static char buf[20];

		char *pos = buf;

		if (x < 0)

			putchar('-'), x = -x;

		do

			*pos++ = x % 10 + '0';

		while (x /= 10);

		while (pos > buf)

			putchar(*--pos);

	}

	inline void write(const double a, const int fix = 9)

	{

		printf("%.*f", fix, a);

	}

	const int B = 17, N = 1 << (B + 1) | 10;

	const double PI = 3.141592653589793238462643383279502884197169399375105820974944L;

	namespace FFT

	{

		int rev[N];

		struct cpx

		{

			double x, y;

			cpx(const double _x = 0.0, const double _y = 0.0)

				: x(_x), y(_y) {}

			cpx operator + (const cpx &b) const

			{

				return cpx(x + b.x, y + b.y);

			}

			cpx operator - (const cpx &b) const

			{

				return cpx(x - b.x, y - b.y);

			}

			cpx operator * (const cpx &b) const

			{

				return cpx(x * b.x - y * b.y, x * b.y + y * b.x);

			}

			cpx conj() const

			{

				return cpx(x, -y);

			}

		}omega[N], inv[N];

		void init(const int lg2)

		{

			for (int i = 0; i < (1 << lg2); i++)

			{

				rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lg2 - 1));

				omega[i] = cpx(cos(2 * PI * i / (1 << lg2)), sin(2 * PI * i / (1 << lg2)));

				inv[i] = omega[i].conj();

			}

		}

		void fft(cpx *a, const int n, const cpx *w)

		{

			for (int i = 0; i < n; i++)

				if (i < rev[i])

					swap(a[i], a[rev[i]]);

			for (int len = 1; len < n; len <<= 1)

				for (int i = 0; i < n; i += (len << 1))

					for (int k = 0; k < len; k++)

					{

						cpx tmp = a[i + k] - w[k * (n / (len << 1))] * a[i + len + k];

						a[i + k] = a[i + k] + w[k * (n / (len << 1))] * a[i + len + k];

						a[i + len + k] = tmp;

					}

		}

		void solve(cpx *a, cpx *b, const int n)

		{

			fft(a, n, omega), fft(b, n, omega);

			for (int i = 0; i < n; i++)

				a[i] = a[i] * b[i];

			fft(a, n, inv);

			for (int i = 0; i < n; i++)

				a[i].x /= n;

		}

	}

	using namespace FFT;

	int n;

	double q[N];

	cpx f[N], g[N], revf[N];

	int work()

	{

		read(n);

		for (int i = 0; i < n; i++)

		{

			read(q[i]);

			f[i] = revf[n - i - 1] = q[i];

		}

		int m = n << 1, lg2 = 0;

		for (n = 1; n < m; n <<= 1)

			++lg2;

		init(lg2);

		for (int i = 0; i < (m >> 1); i++)

			g[i] = (i ? 1.0 / ((double)i * i) : 0.0);

		solve(f, g, n);

		for (int i = 0; i < (m >> 1); i++)

			g[i] = (i ? 1.0 / ((double)i * i) : 0.0);

		for (int i = (m >> 1); i < n; i++)

			g[i] = 0;

		solve(revf, g, n);

		for (int i = 0; i < (m >> 1); i++)

			write(f[i].x - revf[(m >> 1) - i - 1].x), putchar('\n');

		return ~~(0^_^0);

	}

}

int main()

{

	return zyt::work();

}

巴特西

【BZOJ3527】[ZJOI2014] 力（FFT）

题目：

分析：

代码：

最新文章

热门文章