浅谈FFT
Fast Fourier Transportation(FFT)
·多项式的表达
系数表达
对于一个次数界为n的多项式\(A(x)=\sum_{j=0}^{n-1}{a_jx^j}\)而言,其系数表达是由一个系数组成的向量\(a=(a_0,a_1,...,a_{n-1})\)。
点值表达
一个次数界为n的多项式A(x)的点值表达就是一个由n个点值对所组成的集合
\[
{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})}
\]
使得对k=0,1,...,n-1,所有\(x_k\)各不相同,
\[
y_k=A(x_k)
\]
简单的求点值运算我们可以随意代入n个不相同的数,然后得出点对,时间复杂度\(\Theta(n^2)\)。后面可以看到,如果我们用一点巧妙的取值,可以使时间复杂度优化到\(\Theta(n\lg_2n)\)。
求值计算的逆(从一个多项式的点值表达确定的系数表达形式)称为插值。
·多项式运算
\[ C_j=\sum_{k=0}^{j}{A_kB_{j-k}} \]
上述是多项式的乘法,我们把C称为A和B的卷积(convolution),表示成\(C=A\bigotimes B\)。
FFT的主要思路是首先把A和B转成点值表达,然后得到C的点值表达,再逆着做一遍,得到C的系数表达。
·DFT与FFT
上述做法太慢,我们要用\(\Theta(n^2)\)的时间把每个多项式转成点值表达,然后再用\(\Theta(n^2)\)的时间转回来,明显很慢,还不如暴力。
我们想,可不可以使用某些特殊的数,使得每次可以做一次运算就可以得到多个数的呢?答案是有的:单位根复数根
n次单位复数根是满足\(\omega^n=1\)的复数\(\omega\)。n次单位复数根恰好有n个,对于k=0,1,...,n-1,这些根是\(e^{\pi ik/n}\)。为了解释这个表达式,我们用复数的指数形式来定义:
\[
e^{iu}=cos(u)+isin(u)
\]
也就是给定一个单位圆,上面均匀地分布着n个向量,如图:
·关于n次单位复根
以上图为例我们可以发现,每一个n(这里是8)次单位复根都是一个向量,他们在乘法意义下形成一个群。
引理1:(消去引理)
\[
对于任意整数n\geqslant0,k\geqslant0,以及d>0,
\]
\[ \omega^{dk}_{dk}=\omega^{k}_{n} \]
证明
\[
\omega^{dk}_{dk}=(e^{2\pi i/dn})^{dk}=(e^{2\pi i/n})^k=\omega^{k}_{n}
\]
引理2:(折半引理)
\[
如果n>0为偶数,那么n个n次单位根的平方集合就是n/2个n/2次单位根的集合
\]
证明
\[
(\omega^{k+n/2}_{n})^2=\omega^{2k+n}_n=\omega^{2k}_n\omega^n_n=(\omega^k_n)^2
\]
因此\(\omega^k_n\)与\(\omega^{k+n/2}_n\)平方相等。
引理3:(求和引理)
\[
对于任意整数n\geqslant0和不能被n整除的非负整数k,有
\]
\[ \sum_{j=0}^{n-1}(\omega^k_n)^j=0 \]
证明
\[
\sum_{j=0}^{n-1}(\omega^k_n)^j=\frac{(\omega^k_n)^0(1-(\omega^k_n)^n)}{1-\omega^{k}_{n}}=\frac{(\omega^n_n)^k-1}{\omega^{k}_{n}-1}=\frac{(1)^k-1}{\omega^{k}_{n}-1}=0
\]
因为要求k不能被n整除,而且仅当k被n整除时\[\omega^k_n=1\]成立,同时保证分母不为0。
DFT
回顾一下,我们希望计算次数界为n的多项式
\[
A(x)=\sum_{j=0}^{n-1}{a_jx^j}
\]
在\[\omega^0_n,\omega^1_n,...,\omega^{n-1}_n\]处的值。假设A以系数形式给出,接下来定义结果\(y_k\):
\[
y_k=A(\omega^k_n)=\sum_{j=0}^{n-1}a_j\omega^{kj}_n
\]
向量\(y=(y_0,y_1,...,y_{n-1})\)就是系数向量\(a=(a_0,a_1,...,a_{n-1})\)的离散傅里叶变换(DFT)。我们也记作\(y=DFT_n(a)\)。
FFT
通过使用一种称为快速傅里叶变换(FFT)的方法,利用复根的特殊性质,我们就可以在\[\Theta(n\lg n)\]的时间内计算出\(DFT_n(a)\)。
注意:通篇的n我们假设是2的整数次幂。
FFT利用分治策略,采用A(x)中偶数下标的系数与奇数下标的系数,分别定义两个新的次数界为n/2的多项式
\[ A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1} \]
\[ A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1} \]
可以很容易发现
\[
A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)
\]
所以原问题转化为求两个次数界为n/2的多项式\[A^{[0]}(x)\]和\[A^{[1]}(x)\]在点\[(\omega^0_n)^2,(\omega^1_n)^2,...,(\omega^{n-1}_n)^2\]的取值。
所以我们可以发现在求出\[A^{[0]}(x^2)\]和\[A^{[1]}(x^2)\]以后,可以算出两个复根的结果:
\[
y_k=y^{[0]}_k+\omega^k_ny^{[1]}_k
=A^{[0]}(\omega^{2k}_n)+\omega^k_nA^{[1]}(\omega^{2k}_n)
=A(\omega^k_n)
\]
还有
\[
y_{k+(n/2)}=y^{[0]}_k-\omega^{k}_ny^{[1]}_k
=y^{[0]}_k+\omega^{k+(n/2)}_ny^{[1]}_k
=A^{[0]}(\omega^{2k}_n)+\omega^{k+(n/2)}A^{[1]}(\omega^{2k}_n)
\]
\[
=A^{[0]}(\omega^{2k+n}_n)+\omega^{k+(n/2)}A^{[1]}(\omega^{2k+n}_n)
=A(\omega^{k+(n/2)}_n)
\]
所以就有代码:
void FFT(comp *a,int n,int inv){
if(n==1) return;
int mid=n/2;
for (int i=0;i<mid;++i) c[i]=a[i*2],c[i+mid]=a[i*2+1];
for (int i=0;i<n;++i) a[i]=c[i];
FFT(a,mid,inv);
FFT(a+mid,mid,inv);
comp wn={cos(2.0*pi/n),inv*sin(2.0*pi/n)},w={1,0};
for (int i=0;i<mid;++i,w=w*wn){
c[i]=a[i]+w*a[i+mid];
c[i+mid]=a[i]-w*a[i+mid];
}
for (int i=0;i<n;++i) a[i]=c[i];
}
·在单位复数根的插值
现在我们展示如何在单位复数根处插值来完成多项式乘法方案,使得我们把一个多项式从点值表达转换回系数表达。
我们可以把DFT写成矩阵乘积
\[
\left[
\begin{matrix}
y_0\\
y_1\\
y_2\\
\vdots\\
y_{n-1}
\end{matrix}
\right]=
\left[
\begin{matrix}
1 & 1 & 1 & 1 &\cdots& 1\\
1 & \omega_n & \omega^2_n & \omega^3_n &\cdots& \omega^{n-1}_n\\
1 & \omega^2_n & \omega^4_n & \omega^6_n &\cdots& \omega^{2(n-1)}_n\\
1 & \omega^3_n & \omega^6_n & \omega^9_n &\cdots& \omega^{3(n-1)}_n\\
\vdots & \vdots& \vdots& \vdots& \ddots& \vdots\\
1 & \omega^{n-1}_n & \omega^{2(n-1)}_n & \omega^{3(n-1)}_n &\cdots& \omega^{(n-1)(n-1)}_n
\end{matrix}
\right]
\left[
\begin{matrix}
a_0\\
a_1\\
a_2\\
\vdots\\
a_{n-1}
\end{matrix}
\right]
\]
尴尬的是跑得贼慢:
随便卡卡就爆了....
分治难免递归,递归常数大。
于是,考虑改进。
·蝴蝶变换
盗图一张
可以发现,每个下标的二进制形式反过来就是它们最后在序列中的位置。
于是就有了迭代打法。
void FFT(Moon *a,int inv){
int i,j,len;
for (i=0;i<n;++i)
if(i<R[i])swap(a[i],a[R[i]]);
for (len=2;len<=n;len<<=1){
int half=len/2;
Moon w,wn={cos(Pi/half),inv*sin(Pi/half)};
for (j=0;j<n-i;j+=len,w={1,0}){
for (i=0;i<half;++i,w=w*wn){
Moon q=w*a[j+half+i],Q=a[j+i];
a[j+half+i]=Q-q;
a[j+i]=Q+q;
}
}
}
}
int main(){
for (i=0;i<n;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(p-1));
}
原文地址:https://www.cnblogs.com/Chandery/p/11332777.html
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- python实现将range()函数生成的数字存储在一个列表中
- Pytorch 使用不同版本的cuda的方法步骤
- CVE-2020-14645 Weblogic远程代码执行漏洞分析
- PaGoDo:一款功能强大的被动式Google Dork
- 微软轻量级系统监控工具sysmon原理与实现完全分析
- 用Windows电脑训练深度学习模型?超详细配置教程来了
- rConfig中的远程代码执行漏洞分析
- CVE-2020-9964:iOS中的信息泄露漏洞分析
- ReconSpider:一款功能强大的高级OSINT框架
- Python 3.9来了!这十个新特性值得关注
- IRFuzz:一款基于YARA规则的文档文件扫描工具
- 内网渗透测试研究:从NTDS.dit获取域散列值
- 腾讯云大禹高防IP之客户端获取真实IP
- 终极解密输入网址按回车到底发生了什么
- Kafka核心原理的秘密,藏在这 17 张图中