快速傅里叶变换（Fast Fourier Transformation）

最新推荐文章于 2024-09-29 15:06:28 发布

Zerg Wang

最新推荐文章于 2024-09-29 15:06:28 发布

阅读量1.5k

点赞数 1

分类专栏： Algorithm

转载请注明出处，谢谢合作！

本文链接：https://blog.csdn.net/Zerg_Wang/article/details/84992351

版权

Algorithm 专栏收录该内容

12 篇文章 0 订阅

订阅专栏

1.问题引入

当两个多项式f(x)、g(x)相乘时，最直接的解决方法是：枚举f(x)的每一项，与g(x)的每一项分别相乘，再合并同类项，时间复杂度为O(n^2)。实际上，我们可以使用快速傅里叶变换，将时间复杂度优化至O(n*log2n)。

2.背景知识

1.点值表示法

对于一个有n项的多项式，既可以使用n个系数来表示它（这n个系数可表示唯一多项式）：f(x)={a0,a1,a2,...,an−1}

此外，我们也可以使用点值表示：即将n个x的值代入该多项式得到n个f(x)，然后用n个点(x,f(x))也可表示唯一多项式，为什么呢？

如果我们将这n个系数视为未知数，将这n个点的代入f(x)，可以得到n条有n个未知数的方程，那么这些“未知数”也就确定了，“未知数”确定了说明系数确定了，即可表示唯一多项式。

现在我们取n个x，分别代入f(x)和g(x)中，将它们转为点值表示，会得到n对点。再将这n对点的纵坐标（即y值）对应相乘，即可在O(n)的时间内得到f(x)*g(x)这个多项式的点值表达式，最后再转换成系数表达式即可得出结果。

2.复数

复数相乘时，模长相乘，幅角相加。（幅角即复数对应的向量与坐标轴间的夹角）

为什么呢？，根据欧拉公式：

我们可以将复数a+bi先转成极坐标形式：，其中 $\LARGE \mathit{r }$ 为模长， $\LARGE \mathit{\Theta }$ 为幅角。

再转成指数形式：

显而易见，相乘之后，系数相乘（即模长相乘），指数相加（即幅角相加）。

这里讲复数，是因为代入多项式的那n个x不是实数，而是复数。如果将n个实数代入有n项的多项式进行计算，光是转换成点值表达式就要O(n)，而利用复数的性质，我们可以在O(log2n)内转换完毕。

3.离散傅里叶变换

Discrete Fourier Transform，简称DFT。

我们取的这n个复数，其实是在复平面上的一个半径为1的单位圆上取的：

假如要取8个点，则如图，在该单位圆先取第一个点(1,0)，然后逆时针开始等距离取剩余7个点。

以这种方法取出的点有以下特点：（设 $\LARGE W_{n}^{k}$ 为n个点中的第k个）

性质1： $\LARGE W_n^k=(W_n^1)^k$ （由模长相乘，幅角相加的特性可得）

性质2： $\LARGE W_n^k=W_{2n}^{2k}$

性质3： $\LARGE W_n^k=-W_{n}^{k+\tfrac{n}{2}}$ （一个点转半圈（k+n/2）得到的点与原来的点关于坐标原点对称）

性质4： $\LARGE W_{n}^{k}$ 的倒数与其共轭复数相等

注意：这里不要将复平面上的点与点值表示法中的点混淆。这里的点代表一个复数，而点值表示法中的点代表的是两个复数。（x和y，因为x为复数，所以x也是复平面上的点）

4.快速傅里叶变换

Fast Fourier Transformation，简称FFT。

其实DFT的时间复杂度还是O(n^2)，然而，借助计算机，我们可以采用分治的方式将时间复杂度降至O(n*log2n)。

首先对f(x)的奇数项和偶数项进行分离：（假设n为偶数）

$\LARGE \mathit{f(x)=a_{0}+a_{1}x+a_{2}x^{2}+...+a_{n-1}x^{n-1}}$

$\LARGE =a_{0}+a_{2}x^{2}+...+a_{n-2}x^{n-2}+a_{1}x+a_{3}x^{3}+...+$ $\LARGE a_{n-1}x^{n-1}$

$\LARGE =a_{0}+a_{2}x^{2}+...+a_{n-2}x^{n-2}+x(a_{1}+a_{3}x^{2}+...+$ $\LARGE a_{n-1}x^{n-2})$

$\LARGE =f_{0}(x^2)+xf_{1}(x^2)$

代入 $\LARGE W_{n}^{k}$ ：

$\LARGE \mathit{f(W_{n}^{k})=f_{0}(W_{n}^{2k})+W_{n}^{k}f_{1}(W_{n}^{2k})}$

$\LARGE =f_{0}(W_{n/2}}^{k})+W_{n}^{k}f_{1}(W_{n/2}^{k})}$

代入 $\LARGE W_{n}^{k+n/2}$ ：

$\LARGE \mathit{f(W_{n}^{k+n/2})=f_{0}(W_{n}^{2k+n})+W_{n}^{k+n/2}f_{1}(W_{n}^{2k+n})}$

$\LARGE =f_{0}(W_{n}^{2k}W_{n}^{n})+W_{n}^{k+n/2}f_{1}(W_{n}^{2k}W_{n}^{n})}$

$\LARGE =f_{0}(W_{n}^{2k})+W_{n}^{k+n/2}}f_{1}(W_{n}^{2k})}$

$\LARGE =f_{0}(W_{n/2}^{k})-W_{n}^{k}}f_{1}(W_{n/2}^{k})}$

注意：这两条代入式在下面的某两行代码中会直接出现，简直一模一样！

只要得出 $\LARGE f_{0}, f_{1}$ 的值，我们就可以求f(x)，然而 $\LARGE f_{0}, f_{1}$ 也可以继续二分，其二分后的子函数也可继续二分的话……

所以，即使原多项式的项数不为2的次数，一般也会补到2的次数。（多出来的项当作系数为0的项处理即可）

5.离散傅里叶逆变换

Inverse Discrete Fourier Transform，简称IDFT。

现在可以把普通多项式转换为点值表示，那如何从点值表示转为系数表示？

将f离散傅里叶变换后的点值的纵坐标 $\LARGE f(W_{n}^{k})$ 作为另一个多项式 $\LARGE h(x)$ 的系数，然后将以上n个复数的倒数 $\LARGE W_{n}^{-k}$ 代入 $\LARGE h(x)$ 并求其傅里叶变换后其点值的纵坐标，有：（注意：以下的i为实数）

$\LARGE h(W_{n}^{-k})=\sum_{i=0}^{n-1}f(W_{n}^{i})(W_{n}^{-k})^{i}$

$\LARGE =\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_{j}(W_{n}^{i})^{j})(W_{n}^{-k})^{i}$

$\LARGE =\sum_{j=0}^{n-1}a_{j}(\sum_{i=0}^{n-1}(W_{n}^{i})^{j}(W_{n}^{-k})^{i})$ （注意：这里的两个求和函数互换了位置）

$\LARGE =\sum_{j=0}^{n-1}a_{j}(\sum_{i=0}^{n-1}(W_{n}^{j-k})^{i})$

$\LARGE j=k, \sum_{i=0}^{n-1}(W_n^{j-k})^i=n$

$\LARGE j\neq k, \sum_{i=0}^{n-1}(W_n^{j-k})^i=0$ （由等比数列求和可得）

综上， $\LARGE h(W_n^{-k})=na_k$

也就是说，对于点值表示的多项式，我们再求一次它的DFT，不过代入的值是原值的倒数（其实也是原值的共轭复数），然后结果除以n得到结果。

这个过程同样可以用分治来做，此过程也可叫IFFT。

3.代码

#include <cmath>
#include <complex>
#include <stdio.h>
#include <algorithm>

#define PI 3.14159265359
#define MAXN 1100000
using namespace std;

complex <double> a[MAXN],b[MAXN];
int n,m,t,ans[MAXN];
	
void FFT(complex <double> *a,int len,int mode);  //mode表示是FFT还是IFFT，-2表示IFFT
int main()
{
	while (~scanf("%d%d",&n,&m))
	{	
		memset(a,0,sizeof(a));
		memset(b,0,sizeof(b));
		for (int i=0;i<n;i++) 
		{
			scanf("%d",&t);
			a[i].real()=t*1.0;
		}
		for (int i=0;i<m;i++) 
		{
			scanf("%d",&t);
			b[i].real()=t*1.0;
		}
		int len=1;
		while (n+m-1>len) len<<=1;  //两个多项式相乘的结果最多可能有n+m-1项，这里求出能容纳这么多项的最小的一个2的次数
		
		FFT(a,len,2);   
		FFT(b,len,2);  
		
		for (int i=0;i<len;i++) a[i]*=b[i];    //直接逐项相乘，简单粗暴
		 
		FFT(a,len,-2); 
		
	    for (int i=0;i<len;i++)	a[i]/=len;     //IFFT后要逐项除以len
	    
		for(int i=0;i<len;i++) ans[i]=(int)(real(a[i])+0.5);    //四舍五入
		while (ans[len-1]==0) len--;                           //去掉前导零
    	for (int i=0;i<len;i++) printf("%d ",ans[i]);
    	printf("\n");
    }
	return 0;
}

void FFT(complex <double> *a,int len,int mode)
{
	if (len==1) return;
	complex <double> leftArray[len/2];
	complex <double> rightArray[len/2];
	
	for (int i=0;i<len/2;i++)
	{
		leftArray[i]=a[i*2];
		rightArray[i]=a[i*2+1];       //奇偶项分离
	}
	
	FFT(leftArray,len/2,mode);
	FFT(rightArray,len/2,mode);         //递归分治
	
	for (int i=0;i<len/2;i++)    //得益于复数的一些性质，我们计算出a[i]，a[i+len/2]也可得出，因此循环一半即可
	{
		complex <double> w(cos(2*PI*i/len),sin(mode*PI*i/len));  //从单位圆上取代入值，当为IFFT模式时mode=-2,得到的w的虚部也正好要为负
		a[i]=leftArray[i]+w*rightArray[i];            
		a[i+len/2]=leftArray[i]-w*rightArray[i];           //这两行代码简直就是上面那两条公式的翻版
	}
}