博客

*快速傅里叶变换是一种被广泛运用于各个领域算法，尤其是通信领域。它的主要用途是计算两个函数的卷积。在这里我们讨论的是用于多项式乘法的快速傅里叶变换。

一个\(n\)次多项式函数\(f\)有两种表达方式：点值表达法和系数表达法。
系数表达法：
系数表达法是一个多项式函数最直接的表达法。
对于一个多项式函数，它总是可以表达为这样一个形式：
$$f(x)=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+…a_{n}x^{n}$$
点值表达法：
对于一个多项式函数，我们可以得到它的一个取值集合：
$$\{(x_{0},f(x_{0})),(x_{1},f(x_{1})),…(x_{n-1},f(x_{n-1}))\}$$
这个取值集合本身也代表了这个多项式函数。因为，当我们选取好一个值以后，总是可以通过一系列计算利用这个取值集合得到这个值关于这个多项式函数的映射。

暴力计算两个多项式的\(f,g\)的卷积的复杂度根据其表达法的不同有着不同的复杂度。
暴力计算两个系数表达的多项式的卷积的复杂度是\(O(n^2)\)的，具体来说就是将它们像竖式乘法一样相乘。形如：
$$fg=\sum_{i=0}^{n}\sum_{i=0}^{m}a_{i}b_{j}x^{i+j}$$
暴力计算两个取值相同的点值表达的多项式的卷积的复杂度是\(O(n)\)的，具体来说就是将它们的每个点值表达式的每一位的结果相乘。
$$fg={(x_{0},f(x_{0})g(x_{0})),(x_{1},f(x_{1})g(x_{1})),…(x_{n+m},f(x_{n+m})g(x_{n+m})))}$$

在一些应用情况，我们需要用比\(O(n^2)\)更优的复杂度计算两个多项式函数的系数表达法卷积。
很容易可以想到的是将一个多项式函数的系数表达式转换为点值表达式，计算卷积后再转换回去，这样可能就能够以一个可以接受的复杂度计算两个多项式的系数表达式。
复杂度的瓶颈在于点值表达式和系数表达式的相互转换。
快速傅里叶变换(Fast Fourier Transform)泛指是一种通过选取巧妙的取值，使得点值表达式可以快速地和系数表达式相互转换的算法。其中，将系数表达式转换为点值表达式的转换被称为快速离散傅里叶变换，将点值表达式转换为系数表达式的转换被称为快速逆离散傅里叶变换。
为了完成快速傅里叶变换，我们需要了解我们即将选取的一组取值——单位复根，以及它的性质。

方程\(x^n=1\)的所有解构成了\(n\)个单位复根，记作\(\omega_{n}\)，其中的第\(i\)项记作\(\omega_{n}^{i}\)
多项式的单位复根具有一些奇妙的性质：
消去引理：
$$\omega_{dn}^{dk}=e^{\frac{2\pi dki}{dn}}=e^{\frac{2\pi ki}{n}}=\omega_{n}^{k}$$
折半引理：
$$(\omega_{n}^{k+\frac{n}{2}})^2=(e^{\frac{2\pi (k+\frac{n}{2})i}{n}})^2=(e^{\frac{2\pi ki}{n}}*(e^{\pi i}))^2=(-e^{\frac{2\pi ki}{n}})^2=(-\omega_{n}^{k})^2=(\omega_{n}^{k})^2$$
求和引理：
$$\sum_{j=0}^{n-1}(\omega_{n}^{k})^j=\frac{(\omega_{n}^{k})^n-1}{\omega_{n}^{k}-1}=\frac{(e^{\frac{2\pi ki}{n}})^n-1}{\omega_{n}^{k}-1}=\frac{0}{\omega_{n}^{k}-1}=0\ (k>0)$$
折半引理的一个推论：
$$\sum_{k=0}^{n-1}\omega_{n}^{k}=0$$

现在我们考虑如何通过选取单位复根来完成快速地将函数表达在系数表达法和取值表达法之间转换。
我们不由自主地选择分治来解决这个问题。为了实现奇偶分治，我们需要将一个函数的系数表达式扩充到\(2^k-1\)次。
然后我们分别考虑FFT和逆FFT
我们作出了如下定义：
$$f(x)=\sum_{i=0}^{(2^n)-1}a_{i}x^{i}$$
$$f_0(x)=\sum_{i=0}^{2^{(n-1)}}a_{2i}x^{i}$$ $$f_1(x)=\sum_{i=0}^{2^{(n-1)}}a_{2i+1}x^{i}$$
则：
$$f(x)=f_0(x^2)+xf_1(x^2)$$
首先我们考虑如何取值——下面我们证明，当我们已经拥有了一个函数的奇偶分治后的函数的点值表达式，我们可以在线性时间内得到这个函数的点值表达式。
对于某一个取值\(f(\omega_{n}^{k})\ (k\le\frac{n}{2})\)，我们有：
$$f(\omega_{n}^{k})=f_0(\omega_{n}^{2k})+\omega_{n}^{k}f_1(\omega_{n}^{2k})$$
$$f(\omega_{n}^{k})=f_0(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}f_1(\omega_{\frac{n}{2}}^{k})$$
$$f(\omega_{n}^{k+\frac{n}{2}})=f_0(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}f_1(\omega_{n}^{2k+n})$$
$$f(\omega_{n}^{k+\frac{n}{2}})=f_0(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}f_1(\omega_{\frac{n}{2}}^{k})$$
通过这样地递归取值，我们可以在\(O(nlog_2n)\)的时间内将一个多项式从它的系数表达法转化为点值表达法。

然后考虑如何快速地将点值表达法转化为系数表达法。
我们将点值表达法获得的取值作为系数，设出了一个函数\(h\)，它的每一项的系数都是原函数的点值表达式中的取值：
$$f'(x)=\sum_{i=0}^{2^n-1}(a_{i}\sum_{j=0}^{2^n-1}\omega_{n}^{ij}x^{j})$$
据超图灵机所说，我们将\(\omega_{n}^{-k}\)代入这个函数，可以得到如下式子：
$$f'(\omega_{n}^{-k})=\sum_{i=0}^{2^n-1}(a_{i}\sum_{j=0}^{2^n-1}\omega_{n}^{ij}\omega_{n}^{-kj})$$
$$f'(\omega_{n}^{-k})=\sum_{i=0}^{2^n-1}(a_{i}\sum_{j=0}^{2^n-1}\omega_{n}^{j(i-k)})$$
由求和引理可得，当且仅当\(i-k=0\)时，\(\omega_{n}^{j(i-k)}=1\)。其他时候均有\(\omega_{n}^{j(i-k)}=0\)
故有：
$$f'(\omega_{n}^{-k})=a_{k}\sum_{j=0}^{2^n-1}\omega_{n}^{0}=na_{k}$$
故而，我们可以通过将函数\(f’\)通过取每个单位复根的倒数来求得原函数的系数。
而这一过程也可以同样可以通过上述的分治思想来优化。

通过快速离散傅里叶变换将函数从系数表达式转化成点值表达式，然后将点值表达式卷积起来，再将卷积以后的点值表达式通过快速离散傅里叶逆变换转化为系数表达式，这样就在\(O(nlog_2n)\)的时间复杂度内完成了计算多项式函数的系数表达式的卷积。

#include<iostream>
#include<cstdio>
#include<cmath>
#define MID (LEN>>1)
#define PI 3.141592653589793238462

struct Complex{
	double r;
	double i;
	Complex(double rIN=0.0,double iIN=0.0):r(rIN),i(iIN){}
	inline Complex operator+(const Complex &B)const{
		return (Complex){r+B.r,i+B.i};
	}
	inline Complex operator-(const Complex &B)const{
		return (Complex){r-B.r,i-B.i};
	}
	inline Complex operator*(const Complex &B)const{
		return (Complex){r*B.r-i*B.i,r*B.i+i*B.r};
	}
}a[4000005],b[4000005],c[4000005];

inline void FFT(int LEN,Complex *A,Complex *B,int typ){
	if(LEN==1){
		B[0]=A[0];
		return;
	}
//	奇偶项分治 
	for(int i=0;i<MID;++i){
		B[i]=A[i<<1];
		B[i+MID]=A[i<<1|1];
	}
	for(int i=0;i<LEN;++i){
		A[i]=B[i];
	}
	FFT(MID,A,B,typ),FFT(MID,A+MID,B+MID,typ);
	Complex bs(cos((PI*2.0)/LEN),typ*sin((PI*2.0)/LEN)),nw(1.0,0.0);
	for(int i=0;i<MID;++i){
		A[i]=B[i]+nw*B[i+MID];
		A[i+MID]=B[i]-nw*B[i+MID];
		nw=nw*bs;
	}
	for(int i=0;i<LEN;++i){
		B[i]=A[i];
	}
}


int n,m,cnt;
void init(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i){
		scanf("%lf",&a[i].r);
	}
	for(int i=0;i<=m;++i){
		scanf("%lf",&b[i].r);
	}
	cnt=1;
	while(cnt<=n+m){
		cnt<<=1;
	}
	FFT(cnt,a,c,1);
	FFT(cnt,b,c,1);
	for(int i=0;i<cnt;++i){
		a[i]=a[i]*b[i];
	}
	FFT(cnt,a,c,-1);
	for(int i=0;i<=n+m;++i){
		printf("%d ",(int)(a[i].r/cnt+0.5));
	}
}

int main(){
	init();
	return 0;
}