[ZJOI2014]力

嘟嘟嘟


这应该是我的第一篇fft题解吧。


对于这道比较裸的多项式乘法题，主要思路就是推出卷积的式子，然后用fft加速。
卷积我在网上找了半天，只看到一篇博客上写了这个东西。再加上自己的脑补，觉得卷积好像就是这个式子：

[c(i) = sum_{j = 0} ^ {i} a_j *b_{i - j} ]

其中(a, b)是俩多项式，只要符合这个形式，(a, b)就可以用fft加速。（自己意会的，如果不对路过的大佬一定要帮我指出来！）


那么对于这道题，目标就是把题中的式子变成可以卷积的式子。
题意：
求序列(E)，每一项(E_i)满足

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

然后就开始推导了……

[E_i = sum_{j = 0} ^ {i - 1} frac{q_j}{(i - j) ^ 2} - sum_{j = i + 1} ^ {n - 1} frac{q_j}{(j - i) ^ 2} ]

为了和卷积的下标保持统一，这里的下标也从(0)开始。
令(A(i))表示第一个多项式，(B(i))表示第二个多项式，且(f(i) = q_i)，(g(i) = frac{1}{i ^ 2})，规定(g(0) = 0),则

[E_i = A(i) - B(i) ]

[A(i) = sum_{j = 0} ^ {i} f(j) *g(i - j) ]

$B(i) = sum _{j = i} ^ {n - 1} f(j) *g(j - i)$     $= sum _ {j = 0} ^ {n - i - 1} f(i +j) * g(j)$ 如果令$f'(i)$表示$f(i)$翻转后的序列，则 $$B(i) = sum _{j = 0} ^{n - i - 1} f(n - i - j - 1) * g(j)$$ 这样$A(i), B(i)$都符合了卷积的形式，于是fft就可做了。

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<stack>
#include<queue>
#include<complex>
using namespace std;
#define enter puts("") 
#define space putchar(' ')
#define Mem(a, x) memset(a, x, sizeof(a))
#define rg register
typedef long long ll;
typedef double db;
typedef complex<db> cp;
const int INF = 0x3f3f3f3f;
const db eps = 1e-8;
const db PI = acos(-1);
const int maxn = 3e5 + 5;
inline ll read()
{
	ll ans = 0;
	char ch = getchar(), last = ' ';
	while(!isdigit(ch)) {last = ch; ch = getchar();}
	while(isdigit(ch)) {ans = (ans << 1) + (ans << 3) + ch - '0'; ch = getchar();}
	if(last == '-') ans = -ans;
	return ans;
}
inline void write(ll x)
{
	if(x < 0) x = -x, putchar('-');
	if(x >= 10) write(x / 10);
	putchar(x % 10 + '0');
}

int n, len = 1;
cp f1[maxn], f2[maxn], g[maxn], omg[maxn], inv[maxn];

void init()
{
	for(int i = 0; i < len; ++i)
	{
		omg[i] = cp(cos(2 * PI * i / len), sin(2 * PI * i / len));
		inv[i] = conj(omg[i]);
	}
}
void fft(cp* a, cp* omg)
{
	int lim = 0;
	while((1 << lim) < len) lim++;
	for(int i = 0; i < len; ++i)
	{
		int t = 0;
		for(int j = 0; j < lim; ++j) if((i >> j) & 1) t |= (1 << (lim - j - 1));
		if(i < t) swap(a[i], a[t]);
	}
	for(int l = 2; l <= len; l <<= 1)
	{
		int q = l >> 1;
		for(cp* p = a; p != a + len; p += l)	
			for(int i = 0; i < q; ++i)
			{
				cp t = omg[len / l * i] * p[i + q];
				p[i + q] = p[i] - t, p[i] += t;
			}
	}
}

int main()
{
	n = read();
	for(int i = 0; i < n; ++i)
	{
		db q; scanf("%lf", &q);
		f1[i].real(q); f2[n - i - 1].real(q);
		if(i) g[i].real(1.0 / (db)i / (db)i);
	}
	while(len < (n << 1)) len <<= 1;
	init();
	fft(f1, omg); fft(f2, omg); fft(g, omg);
	for(int i = 0; i < len; ++i) f1[i] *= g[i], f2[i] *= g[i];
	fft(f1, inv); fft(f2, inv);	//别忘了再变回系数表示法 
	for(int i = 0; i < n; ++i) 
		printf("%.8lf
", f1[i].real() / len  - f2[n - i - 1].real() / len);
	return 0;
}