嘟嘟嘟
这应该是我的第一篇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)
]
#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;
}