快速傅里叶变换FFT学习笔记

文章目录

• 问题背景

• 前置知识

• • 1.多项式函数的表示法
  • 2.多项式的乘法
  • 3.复数
  • 4.单位根

• 离散傅里叶转换(DFT)

  • 快速离散傅里叶转换(FFT)
  • 快速逆离散傅里叶转换(IFFT)
  • 多项式相乘代码

问题背景

无论是从信息学还是数学的角度来看，多项式都具有重要的应用价值。无论是哪种场景下，在信息学竞赛中频繁涉及利用多项式进行问题求解，并涵盖多项式的各类运算。当两个次数均为n的多项式时，在计算加减法时可以在O(n)时间内完成；相比之下，在处理乘法运算时通常需要O(n^2)的时间复杂度；然而，在所有运算的基础中尤其是快速计算方面，则对乘法提出了更高的要求

前置知识

1.多项式函数的表示法

当多项式 F 被表达成 \sum_{i=0}^na_ix^i 时，则其各a_i被称作该函数的i次项系数或简称为系数。这种表现形式被称为系数表示法。对于某个特定值x_i来说，在多项式中代入它会得到对应的输出值 y_i=F(x_i)。这样就会在平面直角坐标系中确定一个对应点(x_i,y_i)。若给定m+1个互不相同的点(x_0,y_0),(x_1,y_1),\dots,(x_m,y_m), 则可以通过代入这些点来建立m+1个方程组求解该多项式的n+1个未知系数（其中m=n）。由此可见，在这种情况下互不相同的n+1个点足以唯一确定一个n次多项式函数。

2.多项式的乘法

若 F(x)=\sum_{i=0}^na_ix^i，G(x)=\sum_{i=0}^mb_ix^i，那么H(x)=F(x)G(x)=\sum_{i=0}^{n+m}\prod_{j+k=i}a_jb_kx^i这是一个 n+m 次多项式。
如果选取一个点 x_i，则 H(x_i)=H(x)|_{x=x_i}=[F(x)G(x)]|_{x=x_i}=F(x_i)G(x_i)也就是说把两个多项式在 x_i 处的取值相乘就可以得到乘出来的多项式在 x_i 处的取值。这样选取互不相同的 n+m+1 个点就能确定新的多项式了。

可以看出用点值表示法处理问题相较于系数表示法而言更为便捷，并且傅里叶变换即通过将系数表示转换为点值表示来实现这一过程。对于一个n次多项式F(x)，我们需要选择m个特定的点来完成计算。每个点的计算所需时间为O(n)，因此总的时间复杂度为O(nm)。这一复杂度与直接采用系数表示法进行相乘相当，并且由于引入了额外的常数因子而显得相对较高。传统的快速傅里叶变换是对这一算法的一种优化改进。先讲一下复数和单位根的概念

3.复数

在日常生活中，我们主要用到的是实数，而除了实数之外还有另一类数叫做虚数。实数好虚数统称为复数。
复数一般有两种表示方法。若用 c 表示一个复数，那么c=a+bi=re^{i\theta}
事实上，我们可以把一个负数和平面中一个点一一对应，那么两种表示法分别对应这平面直角坐标系和极坐标系的情况。我们称 r=\sqrt{a^2+b^2} 为模长，是一个非负数，称 \theta=\arctan\frac{b}{a} 为辐角，他是一个介于 0 到 2\pi 的数值。
加减法
假设 n=a+bi,m=c+di 为两个任意复数，则 n\pm m=(a+c)\pm(b+d)i
乘法
假设 n=a+bi,m=c+di 为两个任意复数，则 nm=(ac-bd)+(ad+bc)i。如果复数可表示为 n=re^{i\theta},m=\rho e^{i\varphi}，则 nm=(r\rho)e^{i(\theta+\varphi)}。可以看出，两复数相乘，就是模长相乘，辐角相加。
乘方
类似于实数，连乘。

4.单位根

如果复数 x 满足 x^n=1，则称 x 为 n 次单位根，记作 \omega_n，也就是说 \omega_n^n=1。

如图所示，这是一个以 O 圆心，半径为 1 的圆，其中 A,B,\dots,H 为圆上不同的八个点。相邻的两个点之间距离相等。更具前面的知识，我们可以把任一点对应成一个复数 e^{ix\theta}，其中 \theta 为点 A 所对应的辐角，k 为一个 0 到 7 之间的整数，则它的 8 次方为 e^{8ik\theta}=e^{2\pi ik}=1，也就是说，这 8 个复数都为 8 次单位根。类似地，我们可以得到 n 次单位根 \omega_n^k=e^{2\pi ik/n}(k\ 为整数)。
一些性质
1.\omega_n^k=\omega_{n/2}^{k/2}
2.\omega_n^k=\omega_{n}^{k\pm n}
3.\omega_n^{k\pm \frac n 2}=-\omega_n^k
关于这些性质的证明，读者可以根据上面的图来理解。

傅里叶变换(DFT)

该算法将系数表达式转换为基于点值的表达式。具体而言，在将一个次数为n-1的多项式转换为包含n个点的点值表示时，则需将n次单位根逐一代入从而得到所需的结果。亦即该过程可表示为：对于每个k从0到n-1计算对应的函数值...] 其中函数值由a_0, a_1, \dots, a_{n-1}以及\omega_n^k各次幂之和决定。

快速傅里叶变换(FFT)

快速傅里叶变换即为傅里叶变换的一种优化方法

    void FFT(complex *a){
    	complex w,w_n,x,y;  //复数
    	for(int i=0;i<n;i+=1){
    		if(i<rev[i]) swap(a[i],a[rev[i]]);
    	}
    	for(int i=1;i<n;i*=2){  //长度
    		w_n.real=cos(pi/i);
    		w_n.imag=sin(pi/i);  //w_n^1
    		for(int j=0;j<n;j+=i*2){   //数组中每个多项式的起始位置
    			w=1;
    			for(int k=j;k<i+j;k+=1){  //具体位置k-j才是上文的k
    				x=a[k]; y=a[k+i]*w;
    				a[k]=x+y; a[k+i]=x-y;  //(1)(2)式
    				w=w*w_n;  //求出 w_n^(k+1)
    			}
    		}
    	}
    	return;
    }

快速傅里叶逆变换(IFFT)

快速傅里叶逆变换就是把点值表达式转为系数表达式的算法，他可以由傅里叶变换的式子推导得到。考虑 b_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i，把上文傅里叶变换的式子 y_k=\sum_{i=0}^{n-1}a_i(\omega_n^k)^i 代入，得到：
b_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i \\ =\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i \\ =\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i
考虑 c_j=\sum_{i=0}^{n-1}(\omega_n^{j-k})^i。
当 j=k 时c_j=\sum_{i=0}^{n-1}1^i=n
当 j\not=k 时，不妨令 l=j-k\not=0，即 c_j=\sum_{i=0}^{n-1}(\omega_n^l)^i。将 c_j 乘上 \omega_n^l，则 \omega_n^lc_j=\sum_{i=1}^n(\omega_n^l)^i 两式相减，有 c_j=\frac{(\omega_n^l)^n-1}{\omega_n^l-1}=0故 b_k=\sum_{i=0}^{n-1}a_jc_j=na_k因此 a_k=\frac{b_k}{n}，傅里叶逆变换的公式为 a_k=\frac{1}{n}\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i
类似于快速傅里叶变换，我们同样可以得到快速傅里叶逆变换的代码

    void IFFT(complex *a){
    	complex w,w_n,x,y;
    	for(int i=0;i<n;i+=1){
    		if(i<rev[i]) swap(a[i],a[rev[i]]);
    	}
    	for(int i=1;i<n;i*=2){
    		w_n.real=cos(pi/i);
    		w_n.imag=-sin(pi/i);  //-k次方
    		for(int j=0;j<n;j+=i*2){
    			w=1;
    			for(int k=j;k<i+j;k+=1){
    				x=a[k]; y=a[k+i]*w;
    				a[k]=x+y; a[k+i]=x-y;
    				w=w*w_n;
    			}
    		}
    	}
    	for(int i=0;i<n;i+=1) a[i]/=n;  //别忘了除以n
    	return;
    }

多项式乘法代码

基于 FFT 和 IFFT 的实现, 我们便能够完成多项式乘法运算. 具体而言, 在对长度进行枚举时, 时间复杂度为 O(\log_2n); 而在对位置进行枚举时, 则是 O(n) 的复杂度. 因此, 整个算法的时间复杂度仅为 O(n\log_2n).

    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<cmath>
    using namespace std;
    int n,m,len;
    int F,G;
    int rev[2500001],ans[2500001];
    double pi=acos(-1);
    struct complex{
    	double real,imag;
    	void operator=(double x){
    		real=x; imag=0.0;
    		return;
    	}
    	void write(){
    		printf("%lf+%lfi\n",real,imag);
    	}
    }f[2500001],g[2500001];
    complex operator+(complex x,complex y){
    	complex z;
    	z.real=x.real+y.real;
    	z.imag=x.imag+y.imag;
    	return z;
    }
    complex operator-(complex x,complex y){
    	complex z;
    	z.real=x.real-y.real;
    	z.imag=x.imag-y.imag;
    	return z;
    }
    complex operator*(complex x,complex y){
    	complex z;
    	z.real=x.real*y.real-x.imag*y.imag;
    	z.imag=x.real*y.imag+x.imag*y.real;
    	return z;
    }  //以上是一些关于复数的运算
    void FFT(complex *a,int flag){  //FFT/IFFT 模板
    	complex w,w_n,x,y;
    	for(int i=0;i<n;i+=1){
    		if(i<rev[i]) swap(a[i],a[rev[i]]);
    	}
    	for(int i=1;i<n;i*=2){
    		w_n.real=cos(pi/i);
    		w_n.imag=flag*sin(pi/i);
    		for(int j=0;j<n;j+=i*2){
    			w=1;
    			for(int k=j;k<i+j;k+=1){
    				x=a[k]; y=a[k+i]*w;
    				a[k]=x+y; a[k+i]=x-y;
    				w=w*w_n;
    			}
    		}
    	}
    	return;
    }
    int main(){
    	scanf("%d%d",&n,&m);
    	memset(f,0,sizeof(f));
    	memset(g,0,sizeof(g));
    	memset(rev,0,sizeof(rev));
    	for(int i=0;i<=n;i+=1) scanf("%d",&F),f[i]=F;
    	for(int i=0;i<=m;i+=1) scanf("%d",&G),g[i]=G;
    	len=n=n+m; m=1;
    	while((1<<m)<=n) m+=1; n=(1<<m);
    	for(int i=1;i<n;i+=1){
    		rev[i]=((rev[i>>1]>>1)|((i&1)<<(m-1)));
    	}  //这些都是预处理
    	FFT(f,1); FFT(g,1);
    	for(int i=0;i<=n;i+=1) f[i]=g[i]*f[i];  //逐点相乘
    	FFT(f,-1);
    	for(int i=0;i<=len;i+=1){
    		printf("%.0lf ",(f[i].real+0.5)/n);  //系数为整数，保留整数位
    	}
    	printf("\n");
    	return 0;
    }

谢谢观看，记得点赞