学习笔记第三十七节：快速傅里叶变换

最新推荐文章于 2024-12-23 20:56:25 发布

原创最新推荐文章于 2024-12-23 20:56:25 发布 · 1k 阅读

1 ·

CC 4.0 BY-SA版权

学习笔记同时被 2 个专栏收录

61 篇文章

订阅专栏

数论

49 篇文章

订阅专栏

本文深入探讨了FFT算法原理，包括复数运算、单位根性质及循环卷积问题的解决方法，详细讲解了DFT与IDFT过程，提供了高效的递归与循环模拟实现。

正题

解决问题：

for(int i=0;i<n;i++) for(int j=0;j<n;j++) h[i+j]+=f[i]*g[j];

前置知识

定义 $i=\sqrt{-1}$ ，其中i为负数单位。

在这里我们假设z是一个复数，也就是说 $z=a+bi(a,b\in R)$

$Re(z)=a,Lm(z)=b$ ，Re为实部，Lm为虚部。

加减法就是对应相加减：也就是说 $z_1\pm z_2=(a_1\pm a_2)+(b_1\pm b_2)i$

但是乘法有点特殊 $z_1*z_2=(a_1+b_1i)+(a_2+b_2*i)\\=a_1*a_2-b_1*b_2+a_1*b_2i+a_2*b_1i\\=(a_1*a_2-b_1*b_2)+(a_1*b_2+a_2*b_1)i$

共轭 $Z=a-bi$ 。

模长 $\begin{vmatrix} z \end{vmatrix}=\sqrt{a^2+b^2}$

复数除法可以改写成： $z_0/z_1=(z_0*Z_1)/(z_1*Z_1)=(z_0*Z_1)/\ \begin{vmatrix} z_1 \end{vmatrix}^2$

复数 $z=a+bi$ 与复平面上 $(a,b)$ 一一对应，Z与z关于x轴对应。

设 $r=\begin{vmatrix} z \end{vmatrix}$ ， $\varphi$ 为幅角（即向量与x正半轴逆时针旋转的夹角），那么可以得到 $z=r*\cos \varphi +r*\sin \varphi* i$ 。

我们用一种简便的表达方式来表示，就是 $e^{i\varphi}=\cos \varphi +\sin \varphi *i$ 。

我们来讨论一下复数的乘法在复平面上的变化。

首先讨论实部变化：

$\\Re(z_1*z_2)=r_1*\cos \varphi_1*r_2*cos\varphi_2-r_1*\sin \varphi_1*r_2*\sin\varphi_2\\=r_1*r_2*(\cos \varphi_1*\cos \varphi_2-\sin\varphi_1*\sin\varphi_2)\\=r_1*r_2\cos(\varphi_1+\varphi_2)$

其中第二个等号到第三个等号就是两角和公式。

在讨论虚部：

$\\Lm(z_1*z_2)=r_1*\cos\varphi_1*r_2*\sin\varphi_2+r_2*\cos\varphi_2*r_1*\sin\varphi_1\\=r_1*r_2*sin(\varphi_1+\varphi_2)$

也是两角和公式。

那么，乘法就可以写成 $(r_1*e^{i\varphi_1})*(r_2*e^{i\varphi_2})=r_1*r_2*e^{i(\varphi_1+\varphi_2)}$ 。

就是说，两复数相乘的结果，模长等于两者模长之积，幅角为两者之和。

拓展： $z^n=r^n*e^{i(n\varphi)}$

设 $x=r*(\cos\varphi+\sin\varphi*i)$ ，若模长为n次单位根，那么：

$\begin{Bmatrix} r^n=1\\cos(n\varphi)=1 \\ sin(n\varphi)=0 \end{Bmatrix}$

那么我们可以发现，当 $r=1,\varphi=\frac{2\pi k}{n}(k\in[0,n-1])$ ， $\pi$ 为180度，恰好是 $x^n=1$ 的n个不同的解。

在复平面上也很好理解，就是一个单位圆，从原点开始平均分成n份，每一份之间的切割线就是一个解，因为那些切割转n次一定会转到 $(1,0)$

那么 $x=r*\cos\frac{2*\pi*k}{n}+\sin\frac{2*\pi*k}{n} * i,(k\in [0,n-1])$ 。

我们设 $\omega_n=\cos\frac{2*\pi}{n}+\sin\frac{2*\pi}{n}i$ 。

那么n个解就是 $\omega_n^k(k\in[0,n-1])$ 。

现在我们把这个 $\omega_n$ 叫做n次单位根。

单位根的几个性质

首先 $\omega_n^m=\omega_{n/2}^{m/2}$ ，因为相当于两个分数上下都同时除以了2（那首先要保证n,m都是偶数。

其次 $\omega_n^2=\omega_{n/2}$ 这也相当于分数上下同时除以2.

再来 $\omega_n^m=-\omega_n^{m+n/2}$ ，因为当n为偶数的时候 $\omega_n^{n/2}=-1$ ，那么提出来就变成这个样子。

接着是一个最重要的性质：

$\sum_{k=0}^{n-1}\omega_n^{kt}=\begin{Bmatrix} n,n\mid t\\ 0,n\not\ \mid t \end{Bmatrix}$

当 $w_n\neq 1$ 时，很明显是一个等比数列求和，就等于 $\frac{\omega_n^{tn}-1}{\omega_n-1}=0$ 。

否则 $\sum_{k=0}^{n-1}\omega_n^{kt}=\sum_{k=0}^{n-1}1^{kt}=n$ 。

我们就得到了 $\frac{1}{n}\sum_{k=0}^{n-1}\omega_n^{kt}=\begin{Bmatrix} 1,n\mid t\\ 0,n\not\ \mid t \end{Bmatrix}$ 。

前置知识就到这里了。

我们现在开始讲FFT

我们先解决一个叫做循环卷积的东西，

也就是这样的问题：

for(int i=0;i<n;i++) for(int j=0;j<n;j++) h[(i+j)%n]+=f[i]*g[j];

那么我们可以得到一条公式：

$\\h[z]=\sum_{x=0}^{n-1}\sum_{y=0}^{n-1}f[x]*g[y]*[x+y-z\mid n] \\=\sum_{x=0}^{n-1}\sum_{y=0}^{n-1}f[x]*g[y]*\frac{1}{n}\sum_{k=0}^{n-1}\omega_n^{k*(x+y-z)} \\=\frac{1}{n}\sum_{k=0}^{n-1}\omega_n^{-kz}*(\sum_{x=0}^{n-1}\omega_n^{kx}*f[x])*(\sum_{y=0}^{n-1}\omega_n^{ky}*g[y])$

就是用上面单位根的性质将 $[x+y-z\mid n]$ 替换掉。

然后就得出了最后一个东西.

我们设 $F[k]=\sum_{x=0}^{n-1}\omega_n^{kx}*f[x],G[k]=\sum_{x=0}^{n-1}\omega_n^{kx}*g[x],\\H[k]=F[k]*G[k],h[z]=\frac{1}{n}\sum_{k=0}^{n-1}*\omega_n^{-kz}*H[k]$ 。

那么很明显发现求解 $F,G$ 函数的过程是类似的，我们把它称作DFT。

也发现计算H函数是一个简单的 $O(n)$ 问题，计算h函数与计算F，G函数的区别就在于多了-号和 $\frac{1}{n}$ 。

我们把计算h函数的过程叫做IDFT.

也就是说DFT就是：给定f，求 $F[t]=\sum_{x=0}^{n-1}f[x]*\omega_n^{tx}$ 。

暴力 $O(n^2)$ ，考虑特殊性质。

我们设 $fa=f[0,2,...,n-2],fb=f[1,3,...,n-1]$ 。

那么当n为偶数的时候，两个数组的大小都是 $\frac{n}{2}$ 的。

那么他们的F函数又是什么呢？

$\\Fa[t]=\sum_{a=0}^{n/2-1}fa[a]*\omega_{n/2}^{at} \\Fb[t]=\sum_{b=0}^{n/2-1}fb[b]*\omega_{n/2}^{bt}$

我们看一下 $F,Fa,Fb$ 三者的关系。

$\\F[t]=\sum_{a=0}^{n/2-1}f[2a]*\omega_n^{2at}+\sum_{b=0}^{n/2-1}f[2b+1]*\omega_n^{t*(2b+1)} \\=\sum_{a=0}^{n/2-1}fa[a]*\omega_{n/2}^{at}+\omega_n^t\sum_{b=0}^{n/2-1}fb[b]*\omega_{n/2}^{bt} \\=Fa[t]+\omega_n^tFb[t]$

然而这个东西只有当 $t<n/2$ 的情况下才成立，因为 $Fa,Fb$ 在 $t>=n/2$ 时没有意义。

那么 $t>=n/2$ 怎么解决呢？

还是像那样子，只需要将 $t\to t+n/2$ 就可以了。

$\\F[t+n/2]=\sum_{a=0}^{n/2-1}f[2a]*\omega_n^{(t+n/2)*2a}+\sum_{b=0}^{n/2-1}f[2b+1]*\omega_n^{(t+n/2)*(2b+1)} \\=\sum_{a=0}^{n/2-1}fa[a]*\omega_{n/2}^{at}+\sum_{b=0}^{n/2-1}fb[b]*\omega_{n/2}^{bt}*\omega_n^{n/2}*\omega_{n}^t \\=Fa[t]-\omega_n^tFb[t]$

那么只要我们知道了 $Fa,Fb$ ，我们就可以用 $O(n)$ 的时间求出 $F$ 。

那么求解 $Fa,Fb$ 的过程，就相当于不断地递归，那么最后只用递归 $\log n$ 层就可以了，因为每次都是折半的。

在这里我们注意到，每次的n都要是二的倍数，否则很难做。所以，我们令n为大于两多项式最高指数之和的二次幂就可以了。

当只剩下一个数时，明显可以直接返回，因为 $F[0]*\omega_1^0=F[0]$ 。

关于IDFT，只要在求F数组的时候，将 $\omega_n\to\omega_n^{-1}$ 就可以了，相当于变成它的共轭复数，因为一个模长为1的复数，它的共轭复数与它本身相乘恰好等于1，所以在复数意义下互为倒数。

递归显得很简单。

但是如果我们每次都这样给 $Fa,Fb$ 数组赋值递归，时间复杂度明显很大，这时候，就可以用循环来模拟递归。

下面给出模板：

#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;

int n,m;
struct complex{
	double x,y;
	complex operator+(complex a)const {return (complex){x+a.x,y+a.y};}
	complex operator-(complex a)const {return (complex){x-a.x,y-a.y};}
	complex operator*(complex a)const {return (complex){x*a.x-y*a.y,x*a.y+y*a.x};}
}a[3000010],b[3000010];
int lim=1,l=0;
int where[3000010];
double co[3000010],si[3000010];
const double Pi=acos(-1.0)*2.0;
complex wn,w,fa,fb;

void dft(complex *now,int idft){
	for(int i=0;i<lim;i++) if(where[i]>i) swap(now[i],now[where[i]]);
	complex wn,w,a,b;
	for(int l=2;l<=lim;l*=2){
		wn=(complex){cos(Pi/l),idft*sin(Pi/l)};
		for(int i=0;i<lim;i+=l){
			w=(complex){1,0};
			for(int x=i,y=i+l/2;y<i+l;x++,y++,w=w*wn){
				a=now[x],b=now[y]*w;
				now[x]=a+b;
				now[y]=a-b;
			}
		}
	}
}

int main(){
	scanf("%d %d",&n,&m);n++;m++;
	for(int i=0;i<n;i++) scanf("%lf",&a[i].x);
	for(int i=0;i<m;i++) scanf("%lf",&b[i].x);
	while(lim<n+m-1) lim*=2,l++,co[lim]=cos(Pi/lim),si[lim]=sin(Pi/lim);
	for(int i=0;i<lim;i++) where[i]=((where[i>>1]>>1) | ((i&1)<<(l-1)));
	dft(a,1);dft(b,1);
	for(int i=0;i<lim;i++) a[i]=a[i]*b[i];
	dft(a,-1);
	for(int i=0;i<n+m-1;i++) printf("%d ",(int)(a[i].x/lim+0.5));
}

循环模拟关键就在于处理出最后的序列，因为位置变得很散乱。

通过数学归纳法可以发现规律，x所在的位置恰好是x在2进制意义下的翻转。

那么翻转我们可以通过递推来完成，接着每一次循环模拟即可。

关于循环模拟，可以先枚举n的大小，其次就是每一个块，最后是不断地整合，在打的时候还是要注意许多细节。