FFT
首先要说明一个误区,很多人认为FFT只是用来处理多项式乘的,其实FFT是用来实现多项式的系数表示法和点值表示法的快速转换的,所以FFT的用处远不止多项式乘。
FFT的前置知识:点值表示法,复数运算,三角函数。
多项式的系数表示法和点值表示法
系数表示法
[A(x)=\sum_{i=0}^n a_i*x^i ]
点值表示法
不妨将A视为关于x的函数,点值表示法就是在A的图像上取n个点,则该多项式可以被这n个点唯一确定。
[a_i=A(x_i)\ A(x)={(x_1,a_1),(x_2,a_2),\dots,(x_n,a_n)} ]
点值表示法有什么好处呢?
我们知道系数表示法下两多项式相乘是(O(n^2)),但在点值表示法下奇迹出现了:
[A(x)={(x_1,a_1),(x_2,a_2),\dots,(x_n,a_n)} \ B(x)={(x_1,b_1),(x_2,b_2),\dots,(x_n,b_n)} \ A(x)B(x)={(x_1,a_1b_1),(x_2,a_2b_2),\dots,(x_n,a_nb_n)} ]
显然这个是可以O(n)实现的。虽然但是,我们几乎不会在计算中用到点值表示法,但这也给了我们一个解决多项式乘的思路。系数转点值,相乘,点值转系数。
又很显然,我们可以随便取n个数往函数里带,可惜这样又使复杂度回到了(O(n^2))。
于是FFT出现了,FFT使我们可以用(O(n\log n))的复杂度将系数转换成一组特殊的点值,并再把点值转回系数。
复数的计算
简单的理解,复数就是实数加虚数,多少都知道点虚数吧,没错,知道点就够了。
[i=\sqrt{-1}\ z_1=a+bi\ z_2=c+di\ z_1+z_2=a+c+(b+d)i\ z_1-z_2=a-c+(b-d)i\ z_1*z_2=ac-bd+(da+bc)i ]
单位根
记住以下性质(感兴趣可以自己推推,就是基础的三角函数)
[\omega_n^k=cosk\frac{2\pi}{n}+sink\frac{2\pi}{n} \ \omega_n^0=\omega_n^n=1\ \omega_{2n}^{2k}=\omega_n^k ]
好了,现在你已经掌握了所有FFT的前置知识了,自己来推推FFT吧。
正式开始FFT
将(\omega_n^0,\dots,\omega_n^{n-1}) 这n个数带入得到点值表示,于是:
[A(x)=a_0+a_1x+a_2{x^2}+a_3{x^3}+a_4{x^4}+a_5{x^5}+ \dots+a_{n-2}x^{n-2}+a_{n-1}x^{n-1}\ A(x)=(a_0+a_2{x^2}+a_4{x^4}+\dots+a_{n-2}x^{n-2})+(a_1x+a_3{x^3}+a_5{x^5}+ \dots+a_{n-1}x^{n-1})\ A_1(x)=a_0+a_2{x}+a_4{x^2}+\dots+a_{n-2}x^{\frac{n}{2}-1}\ A_2(x)=a_1+a_3{x}+a_5{x^2}+ \dots+a_{n-1}x^{\frac{n}{2}-1}\ A(x)=A_1(x^2)+xA_2(x^2)\ ]
我们将(ωkn(k
[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}(\omega_n^{2k+n})=A_1(\omega_n^{2k}\omega_n^n)-\omega_n^kA_2(\omega_n^{2k}\omega_n^n)=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k}) ]
摘自attak的blog
观察这两个式子
[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\ ]
显然(怎么又是显然)只要求出(A_1(\omega_n^{2k}))和(\omega_n^kA_2(\omega_n^{2k}))就可以得出(A(\omega_n^k))和(A(\omega_n^{k+\frac{n}{2}})),而(A_1(\omega_n^{2k}))和(\omega_n^kA_2(\omega_n^{2k}))又可以进一步递归,正是这一点性质使FFT的复杂度达到了优秀的(O(n\log n))
IFFT
前面提到了,系数转点值,相乘,点值转系数,所以我们还需要能在(O(n\log n))复杂度下完成点值转系数,这就是IFFT,快速傅利叶逆变换。
可以将傅里叶变换和傅立叶逆变换的公式表示写出来(我也不会推逆变换)。
A(x) 表示多项式的系数表示法 B(x)表示多项式的点值表示法
[FFT\ B(x)=\sum_{i=0}^{N-1} A(k)\omega_N^{ik}\ IFFT\ A(x)=\frac 1 N \sum _{i=0} ^{N-1} B(k)\omega_n^{-ik}\ ]
IFFT就是在把FFT的(\omega_n^{ik}改成\omega_n^{-ik}),然后再乘个(\frac 1 N)
考虑(\omega_n^k)的几何意义(高一三角函数)可以得到
[\omega_n^{ik}=a+bi\ \omega_n^{-ik}=a-bi ]
所以IFFT只需要在FFT上做一点改动。
千万别忘了最后乘个(\frac 1 N)
实现
个人认为这是所有算法最难也最重要的部分,然而很多blog都是将这部分一笔带过,所以我决定来详细的讲讲。
递归写法
竟然是递归,那必然有递归极限。每次递归多项式的项数剩下一半,只剩一项时(A(x)=a_1)与(x)无关,所以直接返回就好了。
在求出了(A_1(\omega_n^{2k}))和(\omega_n^kA_2(\omega_n^{2k}))后做两次多项式加减可以得出(A(\omega_n^k))和(A(\omega_n^{k+\frac{n}{2}}))
由于多次的递归和数组新建和赋值使递归写法的常数大的出奇,所以我们需要更好的写法。
重难点来了!
递推做法
首先观察这张图。
我们的递归做法就是从上向下将原数列对半拆,再合并。
我们的合并顺序就是从下到上合并,详细的说先合并(a_0和a_4,a_2和a_6,a1和a_5,a_3和a_7,然后a_0,a_4合并好的整体与a_2和a_6合并好的整体再合并\dots)
观察最上和最下层的数列的二进制表示,发现就是将二进制翻转了。我们有这个神仙操作可以快速的翻转二进制。
rev[i]=(rev[i>>1]>>1)|((i&1)<
于是我们可以(O(n))得到最下层的数列,再依次往上推,即不用递归也不用开大量数组,让代码变得飞快!
void fft(cp *a,int n,int inv)
{
int bit=0;
while((1<>1]>>1)|((i&1)<
#include
using namespace std;
typedef complex cp;
const int N=(1<>n>>m;
for(int i=0;i>a[i];
for(int i=0;i>b[i];
for(lim=1;lim>1]>>1)|((i&1)<
Original: https://www.cnblogs.com/29taorz/p/15712742.html
Author: T_X蒻
Title: 很详细的FFT(快速傅里叶变换)概念与实现
原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/644121/
转载文章受原作者版权保护。转载请注明原作者出处!