数学——快速傅里叶变换(FFT)

Shan xizeng

1. 基础知识

快速傅里叶变换，用来求出两个多项式相乘，如果暴力相乘，时间复杂度为(O(n^2 ))，使用快速傅里叶变换，可以优化到(O(n log n))。

准备知识：

多项式：设(A(x))表示一个n次多项式，则(A(x)=a_0x^0+a_1x^1+cdots+a_{n-1}x^{n-1})

多项式的表示方法：

一是用系数表示法，表示为(sum_{i=0}^{n-1}a_ix^i)，二是点值表示法，表示为对于几个具体的(x)对应的(A(x))的值，最少需要(n)个不同的点就能表示唯一一个(n-1)次多项式。

多项式运算：

加法：

如果使用系数表示法，则将各个系数相加，复杂度为(O(n))；

如果使用点值表示法，则将横坐标相同点的纵坐标相加，复杂度相同。

乘法：

如果使用系数表示法，则设得到的多项式为(sum_{i=0}^nc_ix^i)，其中，(c_i=sum_{j+k=i,0leq j,kleq n}a_jb_k)，很显然，时间复杂度为(O(n^2))；

如果使用点值表示法，则将横坐标相同点的纵坐标相乘，复杂度仍不变，为(O(n))。

向量：

物理、几何学意义：同时具有大小和方向的量。向量相加满足平行四边形定则。

复数：

分为实部与虚部，形如(a+bi),详见。

复数单位根：

将复数(a+bi)中a看做横坐标，b看做纵坐标，则复数单位根为到原点距离为1的点。这些点构成的集合为一个圆，称为单位圆。单位圆的n等分点称为n次单位根，将幅角为正且最小的数设为(omega_n)，则n次单位根分别为(omega_n^2,omega_n^3,dots,omega_n^n)。(omega_n=omega_n^n=1)。根据欧拉公式(e^{ heta i}=cos heta+i sin heta)，(omega_n^k=e^{frac{2pi k i}{n}}=cos k cdot frac{2pi}{n}+isin k cdot frac{2pi}{n})。

易知，(omega_n^{k+frac{n}{2}}=-omega_n^k,omega_n^2=omega_{frac{n}{2}})。(可以画图试试)

2. 快速傅里叶变换

快速傅里叶变换的思想主要是分治。

进行快速傅里叶变换，就是将n个n次单位根分别对多项式求出对应的值。

对于多项式A(x),将其n次单位根带入，则可得到(A(omega_n^k)=sum_{i=0}^na_iomega_n^{ki})

我们将其按照奇偶性分组，得到：(A(omega_n^k)=sum_{i=0}^{n/2-1}a_{2i}omega_n^{2ki}+omega_n^ksum_{i=0}^{n/2-1}a_{2i+1}omega_n^{2ki})

由上面的公式得：(A(omega_n^k)=sum_{i=0}^{n/2-1}a_{2i}omega_{n/2}^{ki}+omega_n^ksum_{i=0}^{n/2-1}a_{2i+1}omega_{n/2}^{ki})

并且

[egin{eqnarray*} A(omega_n^{k+frac{n}{2}}) &=& sum_{i=0}^{frac{n}{2}-1}a_{2i}omega_{frac{n}{2}}^{ki}+omega_n^{k+frac{n}{2}}sum_{i=0}^{frac{n}{2}-1}a_{2i+1}omega_{frac{n}{2}}^{ki} \ &=&sum_{i=0}^{frac{n}{2}-1}a_{2i}omega_{frac{n}{2}}^{ki}-omega_n^ksum_{i=0}^{frac{n}{2}-1}a_{2i+1}omega_{frac{n}{2}}^{ki} end{eqnarray*}
]

这样需要带入的值就减少了一半，时间复杂度为(T(n)=2T(n/2)+O(n)=O(nlog n))。

当然，我们还需要将点值表示转化为系数表示，具体实现就是傅里叶逆变换，就是把原来的傅里叶变换里的(omega_n^k)变成(-omega_n^k)带上然后把系数除以n就好了。证明：我不会。。

实现的快速方法：二进制优化，背板子就好了QAQ

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

const int Maxn=1100000;
const double Pi=3.14159265358979323846;

int n,m,r[Maxn],limit=1,l;

struct complex {
	double x,y;
}a[Maxn],b[Maxn];

complex operator + (complex a,complex b) {
	return (complex) {a.x+b.x,a.y+b.y};
}

complex operator - (complex a,complex b) {
	return (complex) {a.x-b.x,a.y-b.y};
}

complex operator * (complex a,complex b) {
	return (complex) {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}

void fft(complex *a,int type) {
	for(int i=0;i<limit;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int mid=1;mid<limit;mid<<=1) {
		complex Wn=(complex) {cos(Pi/mid),type*sin(Pi/mid)};
		for(int j=0,r=mid<<1;j<limit;j+=r) {
			complex w=(complex) {1,0};
			for(int k=0;k<mid;k++,w=w*Wn) {
				complex x=a[j+k],y=w*a[j+k+mid];
				a[j+k]=x+y;
				a[j+k+mid]=x-y;
			}
		}
	}
}

int main() {
	scanf("%d%d",&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(limit<=n+m) {
		limit<<=1;l++;
	}
	for(int i=0;i<limit;i++) r[i]=r[i>>1]>>1|((i&1)<<l-1);
	fft(a,1),fft(b,1);
	for(int i=0;i<limit;i++) a[i]=a[i]*b[i];
	fft(a,-1);
	for(int i=0;i<=n+m;i++) printf("%d",(int)(a[i].x/limit+0.5));
	return 0;
}