转载请注明原创
似乎,我们集训队对FFT的研究者并没有太多深入,那么,我就来做一做这个东西吧。
涉及名词:
FFT(Fast Fourier Transform): 快速傅里叶变换
DFT(Discrete Fourier Transform): 离散傅里叶变换
FT(Fourier Transform): 连续傅里叶变换
卷积公式:通过两个函数f 和g 生成第三个函数
单位复根:满足的复数
欧拉公式:
欧拉公式推导如下:
令
相消引理:对于任意偶数n > 0,有
快速傅里叶解决的问题类型:1、多项式相乘(大整数相乘) 2、求解偏微分方程 3、数字信号处理(这玩意跟我扯不上关系。。)等等
对于多项式乘法来说,朴素算法是。其表现形式为两个多项式相加。即,,
多项式的表现方式一般为两种:
- 系数表示法
- 点值表示法
系数表示法是我们正常使用的方法,所以他的计算复杂度为,而点值表示法为
点值表示法表示对于 我们采用了n个不同的x值以及他对应的值形成的点对组成。即,则该多项式通过下面的形式就可以唯一的表示出来:
通过点值表达式,我们可以很快地计算出和。而多项式C的次数界是2n的。所以,需要把A和B都扩充到2n的次数级。
在计算的过程中,我们需要先明白两个概念:
- 求值 :系数=>点值,即将多项式形式转化为点值表达式。
- 插值 :点值=>系数,即将点值转化为多项式形式。
而通过的时间复杂度实现多项式乘法的步骤由以下几步完成。
- 使次数界增加一倍,使多项式A和B的次数级增加到。(复杂度为)
- 求值:通过对A和B分别应用阶的FFT计算出和的长度为2n的点值表示。这两个点值表示中包含了两个多项式在次单位根处的值。(复杂度为)
- 点乘:把的值和的值逐点相乘,就可以计算出多项式的点值表示,这个表示中包含了在每个次单位根处的值。(复杂度为)
- 插值:只要对2n个点值对应用一次FFT以计算出其逆DFT,就可以构造出多项式的系数表示。(复杂度为)
DFT
现在对于多项式n
我们取n个值来作为n的点值。
那么,我们定义的y如下
则成向量是系数向量离散傅立叶变换,也写作
接下来,我们就要利用FFT的方式来进行计算了。我们假设n是2的幂。也就是求值的过程。
所以我们可以求出与在点
其中,
为了更好的理解,我们看成下面
我们把式子在转化一下,因为
在转化下
所以,就可以利用二分折半的方法,来缩小规模。这样的话,计算n个点值,就可以通过的时间算出来。
在计算的时候,我们要引入复数型。
在完成了求值了之后,就是点乘。这个也很好求的。直接将两个相乘就可以了。
然后是最后一步,插值。我们把DFT写成。矩阵如下:
所以,接下来,这要乘以一个逆矩阵,就可以变成,公式如下(这个求的过程,我就不列出来,要写好多。。。):
所以,求得a数组等于下式:
我们发现与之前求y数组一样。只是多了一些常量系数。这些可以等算出后面的,在出去这些常量系数。
在编写模式,会发现,FFT的整个过程类似于归并排序,其空间复杂度很大。所以,其常量因子还是蛮大的。
在100个数据,n <= 1000时的时间对比:
1. N^2的复杂度
2. NlogN的复杂度
在10个数据,n <= 10000时的时间对比:
1. N^2的复杂度
2. NlogN的复杂度
在5个数据, n <= 100000时的时间对比:
1. N^2的复杂度
2. NlogN的复杂度
可以发现,在系数比较小的情况,FFT的求解方法还是比较慢的,在n越来越大的时候,就会体现出其优越性。当然由于需要像归并排序,开空间,结果,导致,时间没有体现出那么明显。虽然,我估计我的写法还是有待优化,才开始,我使用了Vector来作传递,发现太慢了,然后,看上交的模板,又换成了指针作为传递的东西,才使效率提升了很多。
想必应该还是有很多优化的地方。
以下为我的FFT代码,参考了其他方式,增加了非递归的方式。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 | /************************************************************************* > File Name: fft.cpp > Author: skt > Mail: sktsxy@gmail.com > Created Time: 2014年08月11日 星期一 18时21分19秒 ************************************************************************/ #include <bits/stdc++.h> using namespace std; #define LL long long #define pb push_back #define mp make_pair #define eps 1e-6 #define PI acos(-1.0) #define cp complex<double> #define MAXN 400005 template <typename T> inline T Max(T a, T b) {return a>b?a:b;} template <typename T> inline T Min(T a, T b) {return a<b?a:b;} typedef vector <int> vi; typedef vector <cp> vcp; int N, NA, NB, a[MAXN], b[MAXN]; LL ans[MAXN]; cp A[MAXN], B[MAXN], Y1[MAXN], Y2[MAXN], C[MAXN], Y[MAXN]; /* * 输入说明:输入 NA, NB。分别表示A数组有NA + 1项, B数组有NB + 1项 * 接下来为 NA + 1 个数表示A数组从第NA项~第0项的系数 * 接下来为 NB + 1 个数表示B数组从第NB项~第0项的系数 */ void input() { for (int i = 0; i <= NA; i ++) { scanf("%d", &a[i]); } for (int i = 0; i <= NB; i ++) { scanf("%d", &b[i]); } NA ++, NB ++; reverse(a, a + NA); reverse(b, b + NB); } /* * 将A和B的系数调整相同 * 将多项式A和B的次数级扩充到2*n * 2*n为2的幂次方 */ void Extend() { N = Max(NA, NB); N = 1 << (int)(ceil(log(2.0 * N) / log(2.0) - eps)); for (int i = NA; i < N; i ++) a[i] = 0; for (int i = NB; i < N; i ++) b[i] = 0; NA = NB = N; for (int i = 0; i < N; i ++) { A[i] = cp(a[i], 0); B[i] = cp(b[i], 0); } } /* * 快速傅立叶 * 时间复杂度: O(NlogN) * One 为 1表示顺时针 * One 为-1表示逆时针 */ void FFT(cp *a, int n, int One, cp *y) { if (n == 1) { y[0] = a[0]; return ; } cp w = cp(1, 0), wn = cp(cos(One * 2 * PI / n), sin(One * 2 * PI / n)); cp *a1 = new cp[n / 2], *a2 = new cp[n / 2]; cp *y1 = new cp[n / 2], *y2 = new cp[n / 2]; for (int i = 0; i < n / 2; i ++) a1[i] = a[2 * i], a2[i] = a[2 * i + 1]; FFT(a1, n / 2, One, y1); FFT(a2, n / 2, One, y2); for (int i = 0; i < n / 2; i ++) { y[i] = y1[i] + w * y2[i]; y[i + n / 2] = y1[i] - w * y2[i]; w *= wn; } delete a1; delete a2; delete y1; delete y2; } /* * 调整数组,这样就可以以非递归的形式实现FFT * (a0 a1 a2 a3 a4 a5 a6 a7) * (a0 a2 a4 a6) (a1 a3 a5 a7) * (a0 a4)(a2 a6) (a1 a5)(a3 a7) * (a0)(a4)(a2)(a6) (a1)(a5)(a3)(a7) */ void change(cp y[], int n) { for (int i = 1, j = n / 2; i < n - 1; i ++) { if (i < j) swap(y[i], y[j]); int k = n / 2; while (j >= k) { j -= k; k /= 2; } if (j < k) j += k; } } /* * 非递归的fft方式, */ void fft(cp a[], int n, int One) { change(a, n); for (int i = 2; i <= n; i <<= 1) { cp wn = cp(cos(One * 2 * PI / i), sin(One * 2 * PI / i)); for (int j = 0; j < n; j += i) { cp w = cp(1, 0); for (int k = j; k < j + i / 2; k ++) { cp u = a[k]; cp v = w * a[k + i / 2]; a[k] = u + v; a[k + i / 2] = u - v; w = w * wn; } } } } /* * 点乘 */ void Dot() { for (int i = 0; i < N; i ++) { C[i] = Y1[i] * Y2[i]; } } void work() { input(); Extend(); FFT(A, N, 1, Y1); FFT(B, N, 1, Y2); Dot(); FFT(C, N, -1, Y); /* * Complex 是double,所以,必须要跟精度转化为int */ for (int i = 0; i < N; i ++) { ans[i] = (LL)ceil(Y[i].real() - eps) / N; } while (N > 0 && !ans[N - 1]) { N --; } for (int i = 0; i < N; i ++) { printf("%lld ", ans[i]); } printf("\n"); } int main() { freopen("fft.in", "r", stdin); freopen("fft.out", "w", stdout); while (scanf("%d %d", &NA, &NB) != EOF) { work(); } return 0; } |
由于本篇文章,写了太多latex指令,所以,为了保证文章的访问速度。做题的代码,我打算放在另一篇文章里。希望大家见谅。我会在后面补上。这篇文章写得一般,希望大神们指教不足。
写的很好啊,原创内容顶一个!
话说你的latex用的是什么插件啊,效果比我的好啊~
内容写得很不错啊!!太赞了!
wp-latex,只是我设置成显示为svg
Pingback: 【转载】快速傅里叶变换FFT | Shangke7788
click 赞!
赞!终于明白是怎么回事儿了!