Featured image of post 浅谈快速傅里叶变换

浅谈快速傅里叶变换

封面来源

Pixiv にこまき-44342920

[toc]

前置芝士

1.卷积(其实没啥用)

https://www.bilibili.com/video/BV1JX4y1K7Dr

正文

1.点值表达式

多项式$A = \sum_{i=0}^{n}{a_ix^i}$,可以转换为另一种表示方法,即将$A$转为函数$f(x)=A(x)$

$$f(x)\Longleftrightarrow{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))….(x_n,f(x_n))}$$

2.多项式乘法

我们知道,两个$n$次多项式相乘的复杂度为$O(n^2)$,完全不能满足某些题的要求,所以考虑优化。那如果用点值表达式来做乘法呢?

$$ A(x)=(x_0,a_0),(x_1,a_1)….(x_n,a_n)\B(x)=(x_0,b_0),(x_1,b_1)….(x_n,b_n)\A(x)B(x)=(x_0,a_0b_0),(x_1,a_1b_1)….(x_n,a_nb_n)$$

复杂度只有$O(n)$!!!

那么总结一下算法的大体思路:

1.将多项式$A,B$转为点值表达式

2.计算多项式乘法

3.把结果转回系数表达式

但是转换的过程仍然是$O(n^2)$,所以还需要优化才行.

3.复数优化

定义

若$\omega^n=1$,则$\omega是$n$次单位根,显然,$n$次单位根有$n$个。

而复数单位根在复平面上模都为$1$,所以这$n$个单位根分布在单位圆上

考虑欧拉公式:

$$ e^{ix}=\cos{x}+i\sin{x} $$

取$x = 2\pi$,则可推出$$ e^{2\pi i}=1=\omega^n$$

此时的单位根称作主次单位根,记作

$$\omega_n=e^{\frac{2\pi i}{n}}$$

在下文中,我们记$\omega_n$为n次单位根。

性质

单位根有三个性质:

1.对$\forall n,k,d \in N$,$\omega_n^k=\omega_{dn}^{dk}$

证明:

$$ \omega_{dn}^{dk}=e^{\frac{2\pi dki}{dn}}=e^{\frac{2\pi ik}{n}}=\omega_{n}^k $$

2.对$\forall n \in {x|x>0且x|2}$,$$ (\omega_n^k)^2=\omega_{\frac{n}2}^k $$

证明:

直接展开

$$(\omega_n^k)^2=\omega_n^{2k}=\omega_{\frac{n}2}^k $$

3.对$\forall n \in {x|x>0且k不整除n}$,$\sum\limits_{j=0}^{n-1}{(\omega_n^k)^j} = 0$

证明:

等比数列求和公式可得

$$\sum\limits_{j=0}^{n-1}{(\omega_n^k)^j} = \frac{(\omega_n^k)-1}{w_n^k-1}=\frac{(1)^k-1}{w_n^k-1}=0$$

FFT

DFT:离散傅立叶变换

对多项式$A(x) =\sum_{i=0}^{n}{a_ix^i}$,取其在$\omega_n^0,\omega_n^1,\omega_n^2,….\omega_n^{n-1}$的值,将其转换为点值表达式$(x_0,\omega_n^0),(x_1,\omega_n^1),…(x_{n-1},\omega_n^{n-1})$.这个过程就是DFT。

假设$n$是$2$的幂,即使它不是$2$的幂,我们也可以通过向高次幂补0来添加。

FFT优化DFT

考虑上文中的多项式$A(x) =a_0+a_1x+a_2x^2\cdots+a_{n-1}x^{n-1}$,重定义为次数为$\frac n 2$的小多项式$A^{[0]}(x)=a_0+a_2x+a_4x^2\cdots+a_{n-2}x^{\frac{n}2-1}$,$A^{[1]}(x)=a_1+a_3x+a_5x^2\cdots+a_{n-1}x^{\frac{n}2-1}$,发现$A(x)=A^{[0]}(x^2)+A^{[1]}(x^2)$,然后就将大问题化为了两个小问题,即求$A^{[0]}(x^2)$和$A^{[1]}(x^2)$,那么,上面那个复杂度为$O(n)$就可以简化为$O(\log n)了。

对于另一边,我们可以通过

$$ A(\omega_n^{k+\frac{n}{2}}) = A^{[0]}(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A^{[1]}(\omega_n^{2k+n}) \Longrightarrow A^{[0]}(\omega_n^{2k}\cdot \omega_n^{n})-\omega_n^{k}A^{[1]}(\omega_n^{2k}\cdot \omega_n^{n})$$

得到

$A(\omega_n^{k+\frac{n}{2}})=A^{[0]}(\omega_n^{2k})-\omega_n^{k}A^{[1]}(\omega_n^{2k})$。

这个操作叫做蝴蝶变换

IDFT:逆离散傅立叶变换

这个其实就是将FFT再搞一遍。

迭代版FFT

但是,由于上面的过程如果写递归的话会有非常大的常数,所以还要进行一些修改。

FFT分治的过程可以简单地描述一下:

发现其实可以反向迭代的。

1.取出儿子节点,算出DFT

2.用这一步的DFT替换之前的

3.直到迭代到根节点为止,否则重复1

代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
inline void Init(int lim){
    FOR(i,0,lim)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
inline void FFT(node *s,int lim,double flag){
    for(register int i=0;i<lim;i++) if(i<r[i]) swap(s[i],s[r[i]]);
    for(register int i=1;i<lim;i<<=1){
        node f(cos(pi/i),flag*sin(pi/i));
        for(register int j=0;j<lim;j+=(i<<1)){
            node t(1,0);
            for(register int k=0;k<i;k++,t=t*f){
                node nx=s[j+k];
                node ny=t*s[i+j+k];
                s[j+k]=nx+ny;
                s[i+j+k]=nx-ny;
            }
        }
    }
}

撒花花~~

★,°:.☆( ̄▽ ̄)/$:.°★

所念皆星河
Built with Hugo
主题 StackJimmy 设计

提供全站CDN服务