快速傅里叶变换
FFT是用来计算离散傅里叶变换(DFT)及其逆变换(IDFT)的快速算法。
两个n次多项式直接相乘所需的时间为O(
n2
),而FFT可以将其复杂度降低为O(
nlogn
)。
令A(x) =
∑n−1j=0ajxj
B(x) =
∑n−1j=0bjxj
C(x) =
∑2n−2j=0cjxj
其中
cj=∑jk=0akbj−k
一个n次多项式可以由n+1个点值确定,例如两点确定一条直线,三点确定一个二次函数。
那么我们通过DFT将由系数表达式确定的A,B分别转换为点值表达式以后,可以得到C的点值表达式以后,通过IDFT可以得到C的系数表达式。
例如:
A的点值表达式{
(x0,y0),(x1,y1)⋅⋅⋅(x2n−1,y2n−1)
}
B的点值表达式{
(x0,y′0),(x1,y′1)⋅⋅⋅(x2n−1,y′2n−1)
}
那么C的点值表达式为{
(x0,y0y′0),(x1,y1y′1)⋅⋅⋅(x2n−1,y2n−1y′2n−1)
}
通过精心的选点,我们可以将系数表达式转为点值表达式的求值过程,与将点值表达式转为系数表达式的插值过程变为O( nlogn )。
实现
折半引理:若n>0且n为偶数,那么n个n次单位复数根的平方的集合就是
n2
个
n2
次单位复根的集合。
我们在定义两个次数界为
n2
的多项式
A[0](x)
和
A[1](x)
:
A[0](x)
=
a0+a2x+a4x2+⋅⋅⋅+an−2xn/2−1
A[1](x)
=
a1+a2x+a3x2+⋅⋅⋅+an−1xn/2−1
其中
A[0](x)
包括所有下标为偶数的系数,
A[1](x)
包括所有下标为奇数的系数,所以有A(x) =
A[0](x2)
+
xA[1](x2)
所以我们将求A(x)的点值转化为求次数界n/2处
A[0](x)
和
A[1](x)
的值。
我们由折半引理可以递归的将以上两个多项式在n/2个n/2次单位复根处的值求出,此时我们需要保证n是2的幂。
然后我们将其递归的过程画出来,以0~7为例,我们可以发现,从左到右,其二进制下标分别为:
000,100,010,110,001,101,011,111 将他们反转得到
000,001,010,011,100,101,110,111
所以我们可以将其转化为迭代的版本!
DFT−1n(y)
的结果为:
aj=1n∑n−1k=0ykω−kjn
所以我们只需要用
ω−1n
代替
ωn
,最后将所有结果除以n即可。
下面给出代码:
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <iostream>
using namespace std;
const double pi = acos(-1.0);
const int MAXN = 300003;
struct comp {
double x, y;
comp(double _x = 0, double _y = 0):x(_x), y(_y) {}
comp operator + (const comp& a) const { return comp(x + a.x, y + a.y); }
comp operator - (const comp& a) const { return comp(x - a.x, y - a.y); }
comp operator * (const comp& a) const { return comp(x * a.x - y * a.y, x * a.y + y * a.x); }
};
int rev[MAXN], T;
comp Sin[MAXN], tmp;
void fft(comp *a) {
for(int i = 1; i < T; i++) if(rev[i] > i) swap(a[rev[i]], a[i]);
for(int i = 2, mid = 1, s = (T >> 1); i <= T; mid = i, i <<= 1, s >>= 1) {
for(int j = 0; j < T; j += i) {
for(int k = j, cur = 0; k < j + mid; k++, cur += s) {
tmp = a[k + mid] * Sin[cur];
a[k + mid] = a[k] - tmp;
a[k] = a[k] + tmp;
}
}
}
}
comp A[MAXN];
int n, m;
inline void init() {
for(T = 1; T <= n + m; T <<= 1);
for(int i = 1; i < T; i++) {
if(i & 1) rev[i] = (rev[i >> 1] >> 1) ^ (T >> 1);
else rev[i] = rev[i >> 1] >> 1;
}
Sin[0] = comp(1, 0); Sin[1] = comp(cos(2 * pi / T), sin(2 * pi / T));
for(int i = 2; i < (T >> 1); i++) Sin[i] = Sin[i - 1] * Sin[1];
}
int main() {
scanf("%d%d", &n, &m);
for(int i = 0; i <= n; i++) scanf("%lf", &A[i].x);
for(int i = 0; i <= m; i++) scanf("%lf", &A[i].y);
init();
fft(A);
for(int i = 0; i < T; i++) A[i] = A[i] * A[i];
for(int i = 0; i < (T >> 1); i++) Sin[i].y = -Sin[i].y;
fft(A);
for(int i = 0; i <= n + m; i++) printf("%d%c", (int)(A[i].y / T / 2 + 0.5), i == n + m ? '\n' : ' ');
return 0;
}
于是我们就可以在O( nlogn )的时间内求得多项式乘法的系数表达式。
对于大家都会的高精度乘法,我们也可以看作多项式的乘法,并用FFT优化。
不过有一个小细节,由于折半引理,所以我们要保证次数为2的次幂,如果不是,我们可以直接添加0使其成为2的次幂,然后计算。
对于带通配符的字符串题,也是可以用FFT做的,例如bzoj4503
题目传送门:bzoj4503
题解传送门:http://blog.youkuaiyun.com/Nickwzk/article/details/56495231
通过这道题,说明只要通过构造适当的式子,使其能通过FFT优化,那就可以在 O(nlogn) 的时间复杂度完成带通配符的字符串匹配。
需要注意的事情:
- 范围需要开到2的次幂以上,不是两倍
- 多次求FFT时需要注意初始化的位置,例如VijosP2000 AxBProblem
- 数据量较大时注意精度问题,不能预处理 ω
- 求IDFT时注意用 ω−1n 代替 ωn
- 最后答案需要除以n(这个n在代码中应该是2 * T)
- 最后结果小数点后要四舍五入