【学习笔记】 - 基础数论：快速傅里叶变换

$\text{FFT}$从入门到回头是岸

快速傅里叶变换，即 $\text{FFT(Fast Fourier Transform)}$，用于$O(n\log n)$计算两个多项式之间的卷积

听起来似乎很难，但其实没有那么难理解

在学习快速傅里叶变换之前，我们需要一些简单的前置知识

前置知识：

• 复数

  复数可以表示为 $a+b\text i$ ，其中 $a$ 称作实部 $b\text i$ 称作虚部

  其中$\text i=\sqrt{-1}$

  比起向量更优秀的是，复数可以进行加减乘除，还可以带入多项式

  复数相乘的规则是 模长相乘，幅角相加

  模长就是这个向量的模长（也可以理解为这个点到原点的距离）

  幅角就是 $x$ 轴正方向逆时针旋转到与这个向量共线所途径的角（也可以理解为从原点出发指向 $x$ 轴正方向的射线逆时针旋转到这个点所经过的角）

  • 复数的实现

    • STL自带

      C++ 的 STL 为我们提供了复数

      • 使用方法

        定义complex<double> x;

        这样直接对其进行加减乘除即可

    • 模拟

      因为 STL 中的复数常数较大，而一般$\text{FFT}$的题都卡常

      所以选择手动模拟

      struct Complex{
          double x=0;
          double y=0;
          Complex operator+(const Complex&t)const{
              return {x+t.x , y+t.y};
          } 
          Complex operator-(const Complex&t)const{
              return {x-t.x,y-t.y};
          }
          Complex operator*(const Complex&t)const{
              return {x*t.x-y*t.y,x*t.y+y*t.x};
          }
      }a[N],b[N];
      

  • $n$次单位根

    快速傅里叶变换需要用到的是一些特殊的复数

    在复平面建立一个单位圆(半径为$1$)，将其进行$n$等分

    这 $n$ 个点的横坐标为实部、纵坐标为虚部可以构成 $n$ 个复数

    从点$(1,0)$开始，逆时针将这n个点从0开始编号，第k个点对应的虚数记作$\omega_{n}^{k}$

    我们称 $\omega_{n}^{1}$ 为 $n$ 次单位根

    • 单位根的性质

      • $①$. $\omega_n^i \ne \omega_n^j$

      • $②$. $\omega_n^k = \cos \frac{2k\pi}{n}+\text i \sin \frac{2k\pi}{n}$

      • $③$. $\omega_n^0=\omega_n^n=1$

      • $④$. $\omega_{2n}^{2k}=\omega_{n}^k$

      • $⑤$. $\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$

• 什么叫多项式乘积？

  就是卷积

  设两个多项式 $A(x)=a_0+a_1x+a_2x^2\cdots$ 和 $B(x)=b_0+b_1x+b_2x^2\cdots$

  其卷积就是两个多项式相乘($A(x)B(x)$)

  我们发现直接做再合并同类项就是$n^2$的时间复杂度

  快速傅里叶变换可以在$O(n \log n)$的复杂度求出两个多项式卷积的值

• 多项式的性质

  一个$n$次多项式，我们可以用$n+1$个不同点唯一确定这个多项式

  这意味着我们可以用$n+1$个点来表示多项式，这被称作点表示法

  • 点表示法有什么好处

    如果我们要求$A(x)$和$B(x)$的乘积$C(x)$

    我们可以用点表示法去求乘积

    如果我们要求$C(x_1)$，我们发现$C(x_1)=A(x_1)B(x_1)$

    那么求出来$C$所用的复杂度为$O(n)$

    这样乘法不再是多项式乘法的瓶颈

正文

虽然乘法不再是多项式乘法的瓶颈，但是我们可以发现把多项式转换成点值表示的朴素算法是$O(n^2)$的，另外直接暴力插值转回多项式也是$O(n^2)$

那么我们只要解决这个问题就可以了

• $\text{FFT}$

  • 系数表示法转化为点表示法

    设 $A(x) = a_0 + a_1x + a_2x^2 +...+a_{n-1}x^{n-1}$

    我们需要让奇偶次数分开

    $$\begin{align} A(x) &=a_0+a_1x+\cdots+a_{n-1}x^{n-1}\\ &=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1}) \end{align}$$

    我们让$A_1(x)=(a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1})$

    $A_2(x)=(a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1})$

    那么

    $$A(x)=A_1(x^2)+xA_2(x^2)$$

    把$\omega_{n}^{k}$带入，大力分讨

    • $k \in [0,\frac{n}{2}-1)$

      带入得到

      $$\begin{align} A(\omega_n^k) &=A_1(\omega_n^{2k})+\omega_{n}^{k}A_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^k)+\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k}) \end{align}$$

    • $k \in [\frac{n}{2},n-1)$

      带入得到

      $$\begin{align} A(\omega_n^{k+\frac{n}{2}}) &=A_1(\omega_n^{2k})+\omega_{n}^{k+\frac{n}{2}}A_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^k)-\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k}) \end{align}$$

    那么求出$A(w_n^{k}) k\in (0,n]$区间的所有值就可以处理了

    因为$A(w_n^{k})=A_1(w_{\frac{n}{2}}^k)A_2(w_{\frac{n}{2}}^k)$

    只需要递归求$\frac{n}{2}$的两个区间即可

    最多会递归$\log n$层，每层都是$O(n)$的

    所以复杂度是$O(n \log n)$

  • 点表示法求出系数表达式

    FFT的逆变换怎么搞？

    已有$n$个点$(\omega_n^{k},A(\omega_n^{k}))$

    我们为了方便将其称之为$(x_k,y_k)$

    $$C_k=\sum_{i=1}^{n-1}y_i(\omega_n^{-k})^i$$

    我们可以把这个当做一个关于$y$的多项式，$B(x)=y_0+y_1x+\cdots+y_{n-1}x^{n-1}$

    $$C_k=B(\omega_{n}^{-k})$$

    我们可以再做一次傅里叶变换来求解系数表示式

    但是$C_k$并不是原多项式的系数，需要再除以一个$n$

    $$\begin{align} C_k=\sum_{i=1}^{n-1}y_i(\omega_n^{-k})^i &=\sum_{i=1}^{n-1}(\sum_{j=0}^{n-1}a_j{\omega_n^{i}}^j){\omega_n^{-k}}^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j-k}))^i\\ &=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^{i}) \end{align}$$

    我们把$\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^{i}$可以看做一个多项式

    $S(x)=1+x+x^2\cdots+x^{n-1}$

    这么这里就相当于$S(\omega_{n}^{k})$

    分两种情况来看

    • $①.$ $k \ne 0$

      $$S(\omega_{n}^{k}) \gets \omega_{n}^{k}S(\omega_{n}^{k})$$

      发现这两个式相等

      所以$(1-\omega_{n}^{k})S(\omega_{n}^{k})=0$

      因为$\omega_{n}^{k}\ne 1$

      所以$S(\omega_{n}^{k})=0$

    • $②.$ $k=0$

      $S(\omega_{n}^k)=S(1)=n$

    因此

    $$\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^{i})=na_k$$

• 实现

  因为递归实现的 FFT 常数巨大，所以我们要用迭代实现

  我们对于要计算的去分为奇数项和偶数项

  分成后的多项式发现可以再分

  分到最后，我们发现，每一位对应的系数的二进制表示是这个位置原本系数对应二进制翻转过来

  比如$a_1(001)$和$a_4(100)$

  还有$a_1(0001)$和$a_7(1000)$

• 代码

  本代码为迭代版FFT，已使用蝴蝶操作，手打复数，常数因子可能较小？

  const double PI=acos(-1);
  struct Complex{
  	double x,y;
  	Complex operator+(const Complex&t)const{
  		return {x+t.x , y+t.y};
  	} 
  	Complex operator-(const Complex&t)const{
  		return {x-t.x,y-t.y};
  	}
  	Complex operator*(const Complex&t)const{
  		return {x*t.x-y*t.y,x*t.y+y*t.x};
  	}
  }a[N],b[N];
  int rev[N],bit,tot;
  int n,m;
  void fft(Complex a[],int inv){
  	for(int i=0;i<tot;i++)
  		if(i<rev[i])
  			swap(a[i],a[rev[i]]);
  	for(int mid=1;mid<tot;mid<<=1){
  		auto w1=Complex({cos(PI/mid),inv*sin(PI/mid)});
  		for(int i=0;i<tot;i+=mid*2){
  			auto wk=Complex({1,0});
  			for(int j=0;j<mid;j++,wk=wk*w1){
  				auto x=a[i+j],y=wk*a[i+j+mid];
  				a[i+j]=x+y,
                  a[i+j+mid]=x-y;
  			}
  		}
  	}
  }
  signed main(){
  	// fire();
  	cin(n,m);
  	for(int i=0;i<=n;i++) cin(a[i].x);
  	for(int i=0;i<=m;i++) cin(b[i].x);
  	while((1<<bit)<n+m+1) bit++;
  	tot=1<<bit;
  	for(int i=0;i<tot;i++)
  		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
  	fft(a,1),fft(b,1);
  	for(int i=0;i<tot;i++)
  		a[i]=a[i]*b[i];
  	fft(a,-1);
  	for(int i=0;i<=n+m;i++)
  		cout((int)(a[i].x/tot+0.5)," ");
  }
  

  解释一下最后的 (int)(a[i].x/tot+0.5)

  在最后取出的复数是包含虚部的，我们只看实部即可

  后面的 +0.5 用于四舍五入