【BZOJ3527】[ZJOI2014] 力（FFT）
题目：
分析：
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();
}
