关于FFT与欧拉公式

快速傅里叶变换

by 云岁月书 $ 2020/7/29 $

updata by 云岁月书 \(2020/8/1\)

预备知识

卷积与循环卷积

定义 \(\Large a,b\) 为两个数列， \(\Large n,m\) 为其长度 。

定义卷积：

\[\Large c_k = \sum\limits_{i+j=k}a_ib_j \]

定义循环卷积（只有 \(n = m\)时存在）：

\[\Large c_k = \sum\limits_{(i+j)\space mod \space n =k}a_ib_j \]

欧拉公式

\[\Large e^{\mathrm{i}\theta} = \cos\theta+\mathrm{i}\sin\theta \]

欧拉公式证明

证明可以不用看，与后文没关系，但是几何意义最好还是看一下。

代数证明：

由麦克劳林公式我们知道：

\[\large e^x=1+\sum\limits_{i=1}^{\infty}\frac{x^i}{i!} \]

\[\large\cos{x}=1+\sum\limits_{i=1}^{\infty}(-1)^{i}\frac{x^{2i}}{(2i)!} \]

\[\large\sin{x}=x-\sum\limits_{i=1}^{\infty}(-1)^{i}\frac{x^{2i+1}}{(2i+1)!} \]

将 \(e^x\) 中的 \(x\) 替换为 \(\mathrm{i}x\) 。

我们有：

\[\large \begin{align}e^{\mathrm{i}x}&=1+\sum\limits_{i=1}^{\infty}\frac{(\mathrm{i}x)^i}{i!}\\&=1+\mathrm{i}x-\frac{x^2}{2!}-\frac{(ix)^3}{3!}\dots\\&=1+\sum\limits_{i=1}^{\infty}(-1)^{i}\frac{x^{2i}}{(2i)!}+\mathrm{i}(x-\sum\limits_{i=1}^{\infty}(-1)^{i}\frac{x^{2i+1}}{(2i+1)!})\\&=\cos\theta+\mathrm{i}\sin\theta\end{align} \]

欧拉公式的几合意义：

首先需要了解负数的两种表达方式，坐标法和幅角法。

坐标法就是虚数的基本定义即：\(z = x + \mathrm{i}y\) ,在复平面上记作 \((x,y)\) 。

我们记复平面 \(\mathrm{xOy}\) 上一复数 \(P(x,y)\) ，也可表示为向量 \(\overrightarrow{OP}\) 。

\(\overrightarrow{OP}\) 的长度记作复数\(z\) 的模或绝对值，由此我们可以记向量 \(\overrightarrow{OP} = (\arcsin{\frac{y}{\sqrt{x^2+y^2}}}，\sqrt{x^2+y^2})\) 。

当点 \(P\) 不是原点，向量 \(\overrightarrow{OP}\) 与 \(x\) 轴正向的夹角称为复数的 \(z\) 辐角，记作 \(Argz\) 。辐角的符号规定为：由正实轴依反时针方向转到 \(\overrightarrow{OP}\) 为正，依顺时针方向转到 \(\overrightarrow{OP}\) 为负。

易证幅角法下两复数相乘为幅角相加，模相乘。

复平面上，\(e^{i\theta}\) 表示沿着单位圆从 \(x\) 轴转动 \(\theta\) 角。

证明如下：

对于 \(e\) 我们在实数范围的定义是：

\[\large e = \lim\limits_{n\rightarrow\infty}(1+\frac{1}{n})^n \]

而其在复数上我们的定义是：

\[\large e^{\mathrm{i}}=\lim\limits_{n\rightarrow\infty}(1+\frac{\mathrm{i}}{n})^n \]

取极限，可以证明 \(e^i = (1\space rad,1)\) 。

即 \(e^{i\theta}，\theta\in[0,2\pi]\) 表示复平面上的单位圆。

单位根

\(x^n -1 = 0\) 的 \(x\) 被称作 \(n\) 次单位根。

一般的在复数集内有 \(n\) 个 \(n\) 次单位根，分别记作 \(w_1,w_2,w_3\cdots,w_{n}\) ，由欧拉公式有：

\[\Large w_k = e^{\frac{2\pi\mathrm{i}}{k}} = \cos\frac{2\pi}{k}+\mathrm{i}\sin\frac{2\pi}{k} \]

由欧拉函数的几何意义可知 \(w_n\) 表示复平面上一个单位圆的 \(n\) 等分，如下图即为 \(8\) 次单位根。

若 \(n,k\) 互质，则称 \(w_k\) 为本原单位根，其另一个意义是指 \(\exists\forall x_l,x_r \in \{w_k^i| i \in [0,n)\},l \neq r\) 时有 \(x_l \ne x_r\) ，即从1开始乘上本原单位根，\(n\) 次后恰好生成全部的 \(n\) 次单位根。

多项式乘法

定义多项式乘法：

\[\Large A(x)B(x) = \sum\limits_{i=0}^n\sum\limits_{j=0}^ma_ib_jx^{i+j} = \sum\limits_{i=0}^{n+m}c_kx^k \]

关于朴素多项式乘法，我们很容易可以得到其时间复杂度 \(\Omicron(n^2)\) 。

而且貌似也没有什么办法优化它。

于是我们我们可以转换一下思路，一个 \(n\) 项多项式可以被 \(n+1\) 个点唯一确定。

于是我们有了点值表示法，即 \((x_1,A(x_1)),(x2,A(x2)\dots(x_n,A(x_n)))\) 。

此时 \(A(x)B(x) = \{(x_1,A(x_1)B(x_1)),(x_2,A(x_2)B(x_2))\dots(x_{Max(n,m)},A(x_{Max(n,m)})B(x_{Max(n,m)}))\}\)。

我们再将它变换成多项式即可。

此时其时间复杂度仍是 \(\Omicron(n^2)\)。

离散傅里叶变换

定义 \(a\) 为一个个数列， \(n\) 为其长度 ，令

\[\Large b_k = \sum\limits_{i=0}^{n-1}a_iw_n^{ki} （1） \]

则称 \(b\) 为 \(a\) 的离散傅里叶变换（ \(DFT\) ）,记作 \(\large b = \mathcal{F}(a)\)

其逆变换（\(IDFT\)）为：

\[\Large a_k = \frac{1}{n} \sum\limits_{i=0}^{n-1}b_iw_n^{-ki} （2） \]

证明如下：

带入式 \((1)\) 到 式\((2)\) 中去有

\[\large \begin{align}\sum\limits_{i=0}^{n-1}b_iw_n^{-ki} &= \sum\limits_{i=0}^{n-1}w_n^{-ki}\sum\limits_{j=0}^{n-1}w_n^{ij}a_j \\&= \sum\limits_{j=0}^{n-1}a_j\sum\limits^{n-1}_{i=0}w_n^{-ki}w_n^{ij}\\&=\sum\limits_{j=0}^{n-1}a_j\sum\limits^{n-1}_{i=0}w_n^{i(j-k)}\end{align} \]

对于 \(\large \sum\limits^{n-1}_{i=0}w_n^{i(j-k)}\) 有：

\[\large j = k,\sum\limits^{n-1}_{i=0}w_n^{i(j-k)} = n \]

\[j \ne k \]

\[\large\begin{align}\sum\limits^{n-1}_{i=0}w_n^{i(j-k)}&=\sum\limits^{n-1}_{i=0}(w_n^{j-k})^i\\&=\frac{(w_n^{j-k})^0(1-(w_n^{j-k})^n)}{1-w_n^{j-k}}\\&=\frac{1-(w_n^n)^{j-k}}{1-w_n^{j-k}}\\&=\frac{1-1}{1-w_n^{j-k}}\\&=0\end{align} \]

\[\Large\therefore \frac{1}{n}\sum\limits^{n-1}_{i=0}b_iw_n^{-ki}=\frac{1}{n}na_k=a_k \]

综上 \(IDFT\) 成立，同理可从 \(IDFT\) 推出 \(DFT\) 。

关于循环卷积与 \(DFT\) 我们有如下定理：

\[\Large\mathcal{F}(a*b)=\mathcal{F}(a)·\mathcal{F}(b) \]

\(\Large ·\) 表示逐点乘积，\(a*b\) 表示 \(a,b\) 的循环卷积。

关于证明：

\[\Large\begin{align}\mathcal{F}(a)·\mathcal{F}(b)&=\sum\limits^{n-1}_{i=0}a_iw^{ki}\sum\limits^{n-1}_{j=1}b_jw^{kj}\\&=\sum\limits^{n-1}_{i=1}\sum\limits^{n-1}_{j=1}a_ib_j(w^{i+j})^k\\&=\sum\limits^{n-1}_{L=0}\sum\limits_{i+j=L}a_ib_j\end{align} \]

回到多项式乘法，我们可以通过 \(DFT\) 求出原式点值表示法，相乘后用 \(IDFT\) 求出答案多项式，其时间复杂度仍是 \(\Omicron(n^2)\) ,但是此时结合本原单位根的性质我们可以优化到 \(\Omicron(n\log n)\) 。

快速傅里叶变换（\(FFT\)）

首先考虑一下本原单位根的性质：

\[\Large (w_{2n}^k)^2 = w_n^k \space\space(3) \]

式 \(3\) 证明如下：

由单位根的几何性质我们知道，\(w_{2n}^k\) 表示一个幅角为 \(\frac{k\pi}{n}\) ，模为 \(1\) 的一个向量，它自己平方等于幅角为 \(\frac{2k\pi}{n}\) ，模也为 \(1\) 的一个向量，即 \(w^k_n\) 。

\[w_{2n}^{n+k}=-w_{2n}^k\space\space(4) \]

式 \(4\) 证明如下：

\(w_{2n}^{n+k}=w_{2n}^n\cdot w_{2n}^k\)

由式 \(3\) :

\((w_{2n}^n)^2=w_n^n = 1\) 即 \(w_{2n}^n = \pm 1\) 。

再结合本原单位根性质我们知道当 \(w_{2n}^{2n}=1\) 故 \(w_{2n}^n=-1\)

综上有 \(w_{2n}^{n+k}=w_{2n}^n\cdot w_{2n}^k=-w_{2n}^k\)

由此令 \(n = 2m\) ,我们将$ A（x） $的项按次数的奇偶分类：

\[\begin{align}A(x) &= \sum\limits_{0\le i < m}a_{2i}x^{2i}+\sum\limits_{0\le i < m}a_{2i+1}x^{2i+1}\\&= \sum\limits_{0\le i < m}a_{2i}x^{2i}+x\sum\limits_{0\le i < m}a_{2i+1}x^{2i}\end{align} \]

定义

\[\begin{align}&A_0(x)=\sum\limits_{0\le i < m}a_{2i}x^{i}\\&A_1(x)=\sum\limits_{0\le i < m}a_{2i+1}x^{i}\end{align} \]

那易有：

\[A(x)=A_0(x^2)+xA_1(x^2) \]

蝴蝶操作

结合单位根性质我们有：

\[\begin{align}A(w_n^k)&=A_0((w_n^k)^2)+w_n^kA_1((w_n^k)^2)\\&=A_0(w_m^k)+w_n^kA_1(w_m^k)\\A(w_n^{m+k})&=A_0((w_n^{m+k})^2)+w_n^{m+k}A_1((w_n^{m+k})^2)\\&=A_0(w_m^k)-w_n^kA_1(w_m^k)\end{align} \]

综上我们可以利用分治或者递归来求 \(DFT\) 了，为了快速进行 \(IDFT\) ,我们可以将 \(FFT\) 中用到的 \(w_n\) 替换为 \(w_n^{-1}\) ,最后在乘 \(\frac{1}{n}\)即可。

具体分治操作：

定义 \(A(i)\) 为多项式 \(A\) 第 \(i\) 项次数。

不难发现每 \(j\) 次拆分序列是按照 \(A(i)\) 二进制下第 \(i-1\) 位划分，若为 \(0\) ,则其在拆分后的左序列，反之在右序列。

当然这种操作过于繁琐，所以我们如过对 \(A(i)\) 进行二进制上的序列翻转，那么上面的要求就可以直接当做比较大小。

如下图：

\[(0001)_2,(0010)_2,(0011)_2,(0100)_2 \Rightarrow (0100)_2,(0010)_2,(0001)_2,(0011)_2 \]

​ $$\Large\Downarrow\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\Large\Uparrow$$

\[(1000)_2,(0100)_2,(1100)_2,(0010)_2 \Rightarrow (0010)_2，(0100)_2，(1000)_2，(1100)_2 \]

由此我们可以将其由递归形式变为递推。

最后一些实现上的小技巧：

\[\Large n|N \Rightarrow w^n=w^{\frac{N}{n}} \]

用来快速预处理所有可能用的本原单位根。

\[(a+b\mathrm{i})^2 = (a^2+b^2)+(2ab)\mathrm{i} \]

将 \(A(x)B(x)\) 变为 \(A(x)^2\) ,由三遍 \(FFT\) ,变成两遍 \(FFT\) 。

代码：

struct Complex
{   
    double x,y;

    Complex(const double x_=0,const double y_=0):x(x_),y(y_){}
	//默认赋值为 0 。
    inline Complex operator + (const Complex& b) const{return Complex(x+b.x,y+b.y);}
    inline Complex operator - (const Complex& b) const{return Complex(x-b.x,y-b.y);}
    inline Complex operator * (const Complex& b) const{return Complex(b.x*x-b.y*y,x*b.y+y*b.x);}
};//自定义的复数结构体。

inline void Swap(Complex& a,Complex& b){Complex c = a;a = b;b = c;}
//一定需要注意的一个点，此处Swap不能有返回值，否则在洛谷上无限TLE，RE。
void FFT(Complex *A,const int *rev,const int Limit,const int Type)
{
	for(reg int i = 1; i <= Limit ; ++i) if(rev[i] > i) Swap(A[i],A[rev[i]]);
	//划分。
	reg Complex Wn,W,x,y;
	
	for(reg int mid = 1; mid < Limit ; mid<<=1)
	{
		Wn = Complex(cos(Pi/mid),Type*sin(Pi/mid));
		
		for(reg int i = 0,Len = mid<<1; i < Limit ; i += Len)
		{
			W = Complex(1,0);
			//蝴蝶操作合并。
			for(reg int j = 0; j < mid ; ++j,W = W * Wn;)
			{
				x = A[i|j];y = W*A[i|j|mid];//此处用二进制运算，因为此处加法不会进位。
				
				A[i|j] = x+y;//等价于 A[i+j] = x+y;
				A[i|j|mid] = x-y;//等价于 A[i+j+mid] = x-y;
			}	
		}
	}
	if(Type == -1) for(reg int i = 0; i <= Limit ; ++i)
    				{
        				A[i].x = A[i].x/Limit;
        				A[i].y = A[i].y/Limit;
   					}//如果是逆变换还要乘上1/n。
}

inline void FFT_Pre(int *rev,const int Limit,const int L)
{
	for(reg int i = 1; i <= Limit; ++i) rev[i] = (rev[i>>1]>>1) | ((i&1)<<L); 
}//预处理翻转出的序号所对应的值。