【BZOJ3527】[ZJOI2014] 力（FFT）

题目：

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('
');
		return ~~(0^_^0);
	}
}
 
int main()
{
	return zyt::work();	
}