多项式算法1：FFT（快速傅里叶变换）

前言

算法简介

快速傅里叶变换 (fast Fourier transform), 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。——百度百科

其实是一个用于快速求解多项式之积的算法。
多项式：形如$$A(x)=\sum _{j=0}^{n-1}{a}_j{x}^j$$的式子，展开来就是$$A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$$
如果我们用$A(x)$，$B(x)$表示两个不同的多项式，现在要求解两个数的乘积，朴素求法是直接相乘，时间复杂度为$\Theta (n^2)$（过程请参考高精乘法），当$n=10^5$时会炸。
不过，当我们换一种方法来表示多项式后，结果会不会好一点呢。考虑用离散的点来表示这个函数，根据高斯消元，不难看出，用至少n个不同的点来表示一个n元函数才可以唯一确定这个n元函数。现假设我用点$$(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)),\cdots ,(x_{n-1},A(x_{n-1}))$$来表示n元多项式$A(x)$，用点$$(x_0,B(x_0)),(x_1,B(x_1)),(x_2,B(x_2)),\cdots ,(x_{n-1},B(x_{n-1}))$$来表示n元多项式$B(x)$。
不难看出$A(x)\ast B(x)$(最高为$n-1$次)可用点$$(x_0,A(x_0)\ast B(x_0)),(x_1,A(x_1)\ast B(x_1)),(x_2,A(x_2)\ast B(x_2)),\cdots ,(x_{n-1},A(x_{n-1})\ast B(x_{n-1}))$$表示。
算法时间复杂度为$\Theta (n)$，已经能达到我们的要求了，那是不是说明我们成功完成了任务了呢。不，那还早呢，我们要的是系数表示的多项式，而不是点值表达式，如果按照$1\to n$的方法来取$x_0-x_{n-1}$的值，用秦九韶算法将系数多项式转为点值，用高斯消元将点值转为系数多项式，那么时间复杂度为$\Theta (n^2)$甚至更高。
那么，我们有没有一个好的算法来优化这一过程呢，这就是我们今天要讨论的问题了，这就是传说中的FFT（快速傅里叶变换）。

前置技能

DFT与IDFT
首先，我们要提一下什么是离散傅里叶变换（DFT），离散傅里叶变换就是将一个系数表达式变为一个点值表达式；而将一个点值表达式变化为一个系数表达式称为离散傅里叶逆变换（IDFT）。FFT就是以巨大常数为代价快速进行DFT和IDFT的算法。
复数
这个算法首先要在选未知数上做文章，如果未知数的乘幂一直都是1，那计算就方便许多。但是这样的数我们只能找到1和-1，如何找到更多满足条件的数呢？
这就有赖于数系的扩充了，在初中阶段，数学老师教育我们，$\sqrt{-1}$这个数是不存在的。但在虚数中，数学家用符号$i$来表示$\sqrt{-1}$。复数就是形如$a\ast i+b$的所有数，当$a=0$时复数就是我们熟知的实数。下面要用到复数，如果还有不懂可以自行搜索，我这里就不再深入探讨复数了。
矩阵
在往下看之前最好要了解矩阵的基本概念，以及单位矩阵，逆矩阵的定义，不懂的请自行上网问度娘。
欧拉公式
$$e^{ix}=\cos x+i\sin x$$

正文

对复数了解一点的人都知道，复数可以表示成平面直角坐标系上的一个点，如下图：

该点可表示为$a+b\ast i$。
在上文中，我们试着用数字1到n带入多项式来表示一个n次多项式，这样计算非常费力，我们需要谨慎选好带入的数。根据经验，乘幂为1或-1时非常好算，考虑1和-1。当n较大时，实数已经不能满足我们的需求了，此时，代入一些复数或许可以满足我们的要求。
以下是用结构体实现复数的代码：

const double pi = acos(-1.0);
struct virt{
   
   
	double r , i;
	virt( double r = 0.0 , double i = 0.0 ) {
   
   
		this->r = r;
		this->i = i;
	}
	virt operator + ( const virt &x ) {
   
   
		return virt( r + x.r , i + x.i );
	}
	virt operator - ( const virt &x ) {
   
   
		return virt( r - x.r , i - x.i );
	}
	virt operator * ( const virt &x ) {
   
   
		return virt( r * x.r - i * x.i , i * x.r + r * x.i );
	}
	//复数除法用不到，就不写了
};

单位复根
鉴于上图用坐标系表示复数，我们是否也可以找到一类复数，它们的乘幂都可能为1或-1？
答案是肯定的，如下图：

在该单位圆上的所有复数都符合条件。
复数的乘法在复平面中表现为辐角相加，模长相乘。利用这一点，我们可以轻易地处理出所有复数的$n$次幂。
复数$\omega$满足${\omega }^n=1$称作$\omega$是$n$次单位根，下图包含了所有的4次单位根(图中圆的半径是1)。

同样的，下图是所有的8次单位根。

由此我们不难发现单位根有如下性质：
$${\omega }_{2n}^{2m}={\omega }_n^m$$$${\omega }_n^m=-{\omega }_n^{m+\frac{n}{2}}$$
FFT的具体过程
FFT就是将系数表示法转化成点值表示法相乘，再由点值表示法转化为系数表示法的过程，第一个过程叫做离散傅里叶变换(DFT)，又称求值，第二个过程叫做离散傅里叶逆变换(IDFT)，又称插值。
DFT详细过程
想要求出一个多项式的点值表示发，需要选出$n$个数分别带入到多项式里面。带入一个数复杂度是$\Theta (n)$，那么总复杂度是$\Theta (n^2)$的，这是无法接受的，但我们可以给多项式代入${\omega }_n^0{\omega }_n^1\dots {\omega }_n^{n-1}$来利用$n$次单位根的性质来加速我们的运算。
设$A_0(X)$为$A(X)$偶次项的和，$A_1(X)$为$A(X)$奇次项的和，则：
$$A_0(X)=a_0{x}^0+a_2{x}^1+\cdots +a_{n-1}{x}^{\frac{n}{2}}$$$$A_1(X)=a_1{x}^0+a_3{x}^1+\cdots +a_{n-2}{x}^{\frac{n}{2}}$$因为:$$A({\omega }_n^m)=a_0{\omega }_n^0+a_1{\omega }_n^m+a_2{\omega }_n^{2m}+a_3{\omega }_n^{3m}+\cdots +a_{n-1}{\omega }^{}$$