/ 信号与系统

快速傅立叶变换

快速傅立叶变换其实不难,编码也不复杂;但是一开始接触的太高深,在阅读过大量有关但无用的资料之后,终于总结了一点自己的心得体会,记于此。

既然BZOJ 2017省选模拟赛Day1 T3用了FFT,那我就来谈一谈FFT吧。

去年寒假在《算法导论》上看到了快速傅立叶变换,没有认真学;去年暑假上了快速傅立叶变换的课,没有实现过;直到今天,我才真正意识到我也能写的出来!

想求两个n1n-1次多项式(选n1n-1次是因为这样总共有nn个系数,比较方便编程计算)的乘积,我们一般是O(n2)O(n^2)的做法。但是有的时候nn一大,就比较棘手。我们需要一个更快的方法。

对于一个多项式,给定一个xx,我们就能求出它的值;由拉格朗日插值法(找规律题神器,建议学习)可知,只要给我nn个点,我就能得到一个经过了所有给定点的n1n-1次多项式函数。对于两个多项式的乘积,在某点的函数值等于在该点的两个函数值的乘积:若C(x)=A(x)B(x)C(x)=A(x)*B(x),则(x,C(x))(x,C(x))可以由(x,A(x))(x,A(x))(x,B(x)(x,B(x)得到。

这给我们提供了一种求多项式乘积的新思路:

A(x),B(x)(x,A(x)),(x,B(x))(x,C(x))C(x)A(x),B(x)\rightarrow(x,A(x)),(x,B(x))\rightarrow(x,C(x))\rightarrow C(x)

第一步叫点值,第三步叫插值。中间求乘积的时间复杂度是O(n)O(n)的,所以只要考虑优化点值和插值。

如果说直接取2n2n个点然后求点值和插值,时间复杂度还是O(n2)O(n^2),而且写起来麻烦得多,常数又大,何苦呢。快速傅立叶变换的好处就是能加快点值和插值。方法?选取一些特殊的点,求点值。

nn个点要选一些比较特殊的点,比如说,单位复数根。

A(x)=a0+a1x+a2x2+...+an1xn1A(x)=a_ 0+a_ 1x+a_ 2x^2+...+a_ {n-1}x^{n-1}

A(x)A(x)ωn0\omega_ n^{0}ωn1\omega_ n^{1}ωn2\omega_ n^{2},……,ωnn1\omega_ n^{n-1}处的函数值。ωn0\omega_ n^{0}ωn1\omega_ n^{1}ωn2\omega_ n^{2},……,ωnn1\omega_ n^{n-1}xn=1x^n=1的不同的nn个根。

A(x)A(x)的奇数项和偶数项拆开来:

A0(x)=a0+a2x+a4x2+...+an2xn2A_ 0(x)=a_ 0+a_ 2x+a_ 4x^2+...+a_ {n-2}x^{\frac{n}{2}}

A1(x)=a1+a3x+a5x2+...+an1xn2A_ 1(x)=a_ 1+a_ 3x+a_ 5x^2+...+a_ {n-1}x^{\frac{n}{2}}

那么

A(x)=A0(x2)+xA1(x2)A(x)=A_ 0(x^2)+xA_ 1(x^2)

再说一遍,我们要求A(x)A(x)ωn0\omega_ n^{0}ωn1\omega_ n^{1}ωn2\omega_ n^{2},……,ωnn1\omega_ n^{n-1}处的函数值。它们在A0(x)A_ 0(x)和在A1(x)A_ 1(x)上的函数值可以递归求解。然后呢?怎么变快了?

0i<n2,(ωni+n2)2=ωn2i+n=ωn2i=ωn2i\forall 0\leq i< \frac{n}{2},(\omega^{i+\frac{n}{2}}_ n)^2=\omega^{2i+n}_ n=\omega_ n^{2i}=\omega^{i}_ {\frac{n}{2}}

就是说,递归的时候,我们只需要求一半的点函数值,合并结果的时候就可以求出所有的点的函数值!

T(n)=2T(n2)+O(n)=O(nlogn)T(n)=2*T(\frac{n}{2})+O(n)=O(n\log n)

点值搞定!

观察点值实际上要求的东西:

[y1y2y3...yn]=[111...11ωn1ωn2...ωnn11ωn2ωn2×2...ωn(n1)2...............1ωnn1ωn2(n1)...ωn(n1)(n1)][a0a1a2...an1]\begin{bmatrix} y_ 1\\ y_ 2\\ y_ 3\\ ...\\ y_ n\end{bmatrix}=\begin{bmatrix} 1 & 1 & 1 & ... & 1\\ 1 & \omega_ n^1 &\omega_ n^2& ... & \omega_ n^{n-1}\\ 1 & \omega_ n^{2} & \omega_ n^{2\times 2} & ... & \omega_ n^{(n-1)\cdot 2}\\ ... & ... & ... & ... & ...\\ 1 & \omega_ n^{n-1} & \omega_ n^{2(n-1)} & ... & \omega_ n^{(n-1)\cdot(n-1)} \end{bmatrix} \begin{bmatrix} a_ 0\\ a_ 1\\ a_ 2\\ ...\\ a_ {n-1} \end{bmatrix}

设中间的变换矩阵为VV。点值就是求Y=VAY=VA,而插值则是点值的逆过程:A=V1YA=V^{-1}YVV的逆矩阵不太好推导,但是最终形式十分简单:

Vij1=ωnijnV^{-1}_ {ij}=\frac{\omega_ n^{-ij}}{n}

那就来验证一下吧。我们只需要说明VV1=IV\cdot V^{-1}=I即可。

P=VV1P=VV^{-1},则

Pij=k=0n1VikVkj1P_ {ij}=\sum_ {k=0}^{n-1}V_ {ik}\cdot V^{-1}_ {kj}

=k=0n1ωnikωnkjn=\sum_ {k=0}^{n-1}\omega_ n^{ik}\frac{\omega_ n^{-kj}}{n}

=1nk=0n1ωnk(ij)=\frac{1}{n}\sum_ {k=0}^{n-1}\omega_ n^{k(i-j)}

i=ji=j时,

Pij=1nk=0n1ωn0=1nn=1P_ {ij}=\frac{1}{n}\sum_ {k=0}^{n-1}\omega_ n^0=\frac{1}{n}\cdot n=1

iji\neq j时,

Pij=ωnijnk=0n1ωnkP_ {ij}=\frac{\omega_ n^{i-j}}{n}\sum_ {k=0}^{n-1}\omega_ n^k

=ωnijnωn0ωnn1ωn11ωn1=\frac{\omega_ n^{i-j}}{n}\cdot\frac{\omega_ n^0-\omega_ n^{n-1}\omega_ n^1}{1-\omega_ n^1}

=ωnijnωn0ωnn1ωn1=\frac{\omega_ n^{i-j}}{n}\cdot\frac{\omega_ n^0-\omega_ n^n}{1-\omega_ n^1}

=0=0

由此可以看出PP事实上就是单位矩阵。

那求点值和求插值有什么区别?仅仅是把ωn1\omega_ n^1换成了ωn1\omega_ n^{-1},算出答案之后再/n/n就可以了。

插值完毕。

实现上的优化

把每一层实际操作的元素递归写出来,发现元素ii的位置是在ii二进制倒过来那里。这个叫bit-reverse操作,这么做之后就不需要用递归实现FFT了,只需要用循环归并答案。

细节很复杂,但是写出来还是很简单的。

typedef complex<double> comp;
inline comp omega(int i,int n){
	double angle=M_PI*2/n*i;
	return comp(cos(angle),sin(angle));
}
inline void fft(comp *a,int d){
	static int rev[N];
	rev[0]=0;
	for(int i=0;i<N;i++){
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(SHIFT-1));
		if(i<rev[i]){
			swap(a[i],a[rev[i]]);
		}
	}
	for(int i=2;i<=N;i<<=1){
		int half=i>>1;
		for(int j=0;j<half;j++){
			comp w=omega(d*j,i),p,q;
			for(int k=j;k<N;k+=i){
				p=a[k],q=a[k+half]*w;
				a[k]=p+q,a[k+half]=p-q;
			}
		}
	}
	if(d==-1){
		for(int i=0;i<N;i++){
			a[i]/=N;
		}
	}
}

拓展

因为现在时间比较紧了,最后这一段也只是为了填坑而已,所以傅立叶变换之后的许多性质还有NTT、FWT那些也没有时间记录了。FFT的学习就先这样告一段落吧。