嘟嘟嘟
这应该是我的第一篇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\n", f1[i].real() / len - f2[n - i - 1].real() / len);
return 0;
}
原文:https://www.cnblogs.com/mrclr/p/10088928.html