FFT、NTT小结 A Summary for FFT、NTT

Ⅰ. 概述

• 多项式乘法，大整数乘法，直接计算复杂度是$O(n^2)$ 的，有一种分治乘法，$Karatsuba$乘法，可以做到$O(n^{log_23})$
• 快速傅里叶变换$(Fast\ Fourier\ Transform,\ FFT)$可以在$O(nlogn)$的时间计算出两个多项式的乘积
• 快速数论变换$(Number\ Theoretic\ Transform,\ NTT)$，即是$FFT$的数论版本

Ⅱ. 多项式

1. 定义

形如$a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}的代数表达式称为多项式$
其中$a_0, a_1, a_2, \cdots, a_{n-1}$称为多项式的系数, $x$是未知数
未知数的最高次数$n-1$称为多项式的次数，次数$+1=n$称为多项式的次数界, 即多项式的次数小于次数界
记多项式$A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}$，即:

$$A(x)= \displaystyle\sum_{i=0}^{n - 1} a_ix^i$$

• 多项式值的计算

  可利用霍纳法则(又名秦九韶算法)在$O(n)$时间复杂度完成
  $A(x_0)=a_0+x_0(a_1+x_0(a_2+\cdots+x_0(a_{n-2}+x_0a_{n-1}))$

2. 系数表示法

将$n-1$次多项式$A(x)$的系数, 看作一个$n$维列向量, 即：

$$\overrightarrow{a}=(a_0, a_1, a_2, \cdots, a_{n-1})^{T}$$

3. 点值表示法

选取$n$个不同的数$x_0,x_1,⋯,x_{n-1}$对多项式$A(x)$进行求值,得到$A(x_0), A(x_1), \cdots, A(x_{n-1})$
那么就称$\{\left(x_i, A(x_i)\right), 0 \leq i < n, i \in \mathbb Z\}$为多项式$A(x)$的点值表示

4. 多项式的加法和乘法

• 加法

  如果$A(x)$和$B(x)$是次数界为$n$的多项式，那么和也是一个次数界为$n$的多项式$C(x)$，也就是说：
  

  $$A(x)=\sum_{i=0}^{n-1}a_ix^i, B(x)=\sum_{i=0}^{n-1}b_ix^i, 则C(x)=\sum_{i=0}^{n-1}c_ix^i, 其中c_i=a_i+b_i$$
  这是系数表达式，可以在$O(n)$时间复杂度完成

• 乘法

  如果$A(x)$和$B(x)$分别是次数界为$n_a$和$n_b$的多项式，那么积是一个次数界为$n_a+n_b−1$的多项式$C(x)$，
  也就是说：
  

  $$C(x)=\sum_{i=0}^{2n-2}{c_ix^i}, 其中c_i=\sum_{j=0}^{i}a_jb_{i-j}$$
  这是系数表达式，朴素算法可以在$O(n^2)$时间复杂度完成
  计算出的系统向量$\overrightarrow{c}$也称为输入向量$\overrightarrow{a}$和$\overrightarrow{b}$的卷积，表示成$\overrightarrow{c}=\overrightarrow{a}⊗\overrightarrow{b}$
  如果知道$A(x)$和$B(x)$的点值表达, 我们可以$O(n)$把它们逐点相乘得到$C(x)$的点值表达

5. 插值——点值表示法的作用

求值计算的逆(从一个多项式的点值表示确定其系数表示)称为插值
插值多项式的唯一性：
对于任意$n$个不同点值对组成的集合，其中所有的$x$都不同，那么存在唯一的次数界为$n$的对应的多项式$A(x)$

• 插值的计算方法
  观察下面这个矩阵表达式，它通过系数表达式计算出了点值表达式
  

  $$\left[ \begin{array} {ccccc} 1 & x_0 & x_0^2 & ... & x_0^{n - 1} \\ 1 & x_1 & x_1^2 & ... & x_1^{n - 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n-1} & x_{n-1}^2 & ... & x_{n-1}^{n - 1} \\ \end{array} \right] \left[ \begin{array} {c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \\ \end{array} \right] = \left[ \begin{array} {c} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \\ \end{array} \right]$$
  将这个矩阵表达式简写为$A\cdot B = C$, 其中$A$为范德蒙矩阵，$B$为系数矩阵，$C$为点值矩阵
  观察发现，如果求出$A$的逆矩阵$A^{-1}$，两边同时左乘$A^{-1}$, 即：
  $$A^{-1}\cdot A\cdot B=A^{-1}\cdot C\Rightarrow I\cdot B=A^{-1}\cdot C\Rightarrow B=A^{-1}\cdot C$$
  这样我们就通过点值矩阵$C$计算出了系数矩阵$B$，即插值
  但是现在问题来了，这个矩阵是否一定有逆矩阵呢，证明:
  范德蒙矩阵$A$的行列式值$|A|=\prod\limits_{0 \ \leq \ j \ < \ i \; \leq \; n - 1 } (x_i - x_j)$
  显然当这$n$个点都互不相同的时候，行列式的值不为零，所以逆矩阵存在 $\blacksquare$

  计算逆矩阵和矩阵乘法的时间复杂度都是$O(n^3)$的，这样是不行的

  $PS：$现在就是$FFT$大展神威的时刻了，通过选定特定的值来$O(nlogn)$使得系数表达式和点值表达式之间相互计算

Ⅲ. 快速傅立叶变换

• 快速傅立叶变换计算过程：
  
  • 加倍次数界 加入$n$个系数为$0$的高次系数，将多项式$A(x)$和$B(x)$的次数界变为$2n$，时间复杂度$O(n)$
  • 求值 应用$DFT$计算出$A(x)$和$B(x)$的点值表达，时间复杂度$O(nlogn)$
  • 逐点相乘 通过逐点相乘计算$C(x)$的点值表达，时间复杂度$O(n)$
  • 插值 应用$IDFT$计算出多项式$C(x)$的系数表达，时间复杂度$O(nlogn)$

1. 单位复数根

• 前面说到的选定特定的值其实就是单位复数根
• $n$次单位复数根 满足$ω_n^n$=1的复数$ω_n$，$n$次单位复数根恰好有$n$个，为$e^{2πki/n}, k=0,1,2,\dots,n-1$
  解释 根据 欧拉公式 $e^{i \theta} = cos \theta + i * sin \theta$ 可知，$e^{2πi}=1$
  所以$ω_n=e^{2πi/n}$是一个$n$次单位复数根，称为主$n$次单位根，其他$n$次单位根都是$ω_n$的幂次

• 观察这些$n$次单位复数根，$\theta \in [0, 2\pi)$ ，且这$n$个点的$\theta$是均匀地分配在这个区间内的
  在复平面中， $e^{\theta_1} * e^{\theta_2}$的意义相当于$e^{\theta_1}$这个复数逆时针旋转了$\theta_2$的角度
  这$n$个单位复数根将复平面均匀地分成了$n$个部分，每部分弧度为$2\pi/n$
  所以这$n$个单位复数根各不相同 $\blacksquare$

• $n$个$n$次单位复数根在乘法意义下形成了一个群
  因为$\omega_n^n=\omega_n^0=1, \omega_n^i\omega_n^j=\omega_n^{i+j}=\omega_n^{(i+j)mod~n},\omega_n^{-1}=\omega_n^{n-1}$(可通过共轭复数的性质推导)

2. Cooley-Tukey算法

• $FFT$最常见的算法是 Cooley-Tukey 算法，它的基本思路在 1965 年由 J. W. Cooley 和 J. W. Tukey 提出的
  它是一个基于分治策略的算法

• 虽然前面说的第一步将次数界加倍, 但是下面讲解过程仍将多项式次数界看作$n$

2.1 $DFT$

将$n$个$n$次单位根$\omega_n^0,\omega_n^1,…,\omega_n^{n-1}$代入$A(x)$将其转换成点值表达：

$$A(\omega_n^k) = \sum_{i=0}^{n-1}a_i\omega_n^{ki} , k = 0, 1, \cdots, n - 1$$
点值向量$\vec y=(A(\omega_n^0), A(\omega_n^1), \cdots, A(\omega_n^{n-1}))$称作系数向量$\vec a=(a_0, a_1, \cdots, a_{n-1})$的离散傅里叶变换$(Discrete\ Fourier\ Transform, DFT)$，也记作 $\vec y=DFT_n(\vec a)$

2.2 计算$DFT$

• 直接计算$DFT_n(\vec a)$还是需要$O(n^2)$的时间，Cooley-Tukey 算法接下来做的事情是将每一项按照指数奇偶分类
• 我们定义两个新的次数界为$n/2$的多项式$A^{[0]}(x)$和$A^{[0]}(x)$：
  $$A^{[0]}(x)=a_0+a_2x+a_4x^2+…+a_{n-2}x^{n/2-1}$$
  $$A^{[1]}(x)=a_1+a_3x+a_5x^2+…+a_{n-2}x^{n/2-1}$$
• 对比$A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1},可以得到：$
  $$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$$
• 求$A(x)$在$\omega_n^0,\omega_n^1,\omega_n^2,…,\omega_n^{n-1}$的值，即$DFT_n(\vec a)$的问题转化为：
  
  • 求次数界为$n/2$的多项式$A^{[0]}(x)$和$A^{[0]}(x)$在点$(\omega_n^0)^2,(\omega_n^1)^2,…,(\omega_n^{n-1})^2$的值
  • 利用$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$综合结果
• 如果直接这样递归下去，系数不断减半，但是代入的单位复根还是$n$个，这样的话复杂度还是$O(n^2)$
• 消去引理 任何整数$n⩾0，k⩾0，d⩾0，\omega_{dn}^{dk}=(e^{2\pi i/dn})^{dk}=(e^{2\pi i/n})^k=\omega_n^k$
• 根据消去引理，可得折半引理，$(\omega_n^k)^2=\omega_n^{2k}=\omega_{n/2}^k$
• 由折半引理我们可以发现, 每次代入的值也减半了，因为$n/2$次单位复根只有$n/2$个
• 这样递归就没有问题了，边界条件是只剩一个系数的时候，$y_0 = a_0 \omega^0_1 = a_0 * 1 = a_0$
• 只需要计算$k = 0, 1, 2, ..., n/2 - 1$的情况，$w_n^k$称为旋转因子，接下来开始介绍怎么计算：
  $$\begin{align*} y_k &= y^{[0]}_k + y^{[1]}_k\\ &= A^{[0]}(\omega^{2k}_n) + \omega^{k}_n A^{[1]}(\omega^{2k}_n) \\ &= A(\omega^k_n) \end{align*}$$

$$\begin{align*} y_{k + n / 2} &= A(\omega^{k + n / 2}_n) \\ &= A^{[0]}(\omega^{2k + n}_n) + \omega^{k + n / 2}_n A^{[1]}(\omega^{2k + n}_n) \\ &= A^{[0]}(\omega^{2k}_n) - \omega^k_n A^{[1]}(\omega^{2k}_n) \\ &= y^{[0]}_k - \omega^k_n y^{[1]}_k \end{align*}$$

• 蝴蝶操作 假设$u= y^{[0]}_k$，$t= y^{[1]}_k$，且当前旋转因子为$w_n^k$
  我们只需要算出$u$和$t$，那么$y_k = u + \omega^k_n t$，$y_{k+n/2} = u - \omega^k_n t$

• 这个递归计算的时间复杂度，$T(n)=2T(n/2)+O(n)=O(nlogn)$

• 伪代码：

  function RECURSIVE-FFT(a)
  n = a.length
  if n == 1 return a
  wn = e^(2*pi*i/n)
  w = 1
  a0 = (a_0,a_2,...,a_n-2)
  a1 = (a_1,a_3,...,a_n-1)
  y0 = RECURSIVE-FFT(a0)
  y1 = RECURSIVE-FFT(a1)
  for k = 0 to n/2-1
          y[k] = y0[k] + w * y1[k]
          y[k+n/2] = y0[k] - w * y1[k]
          w = w * wn
  

2.3 $IDFT$

• 系数向量$\vec a=(a_0, a_1, \cdots, a_{n-1})$称作点值向量$\vec y=(A(\omega_n^0), A(\omega_n^1), \cdots, A(\omega_n^{n-1}))$的离散傅里叶逆变换$(Inverse\ Discrete\ Fourier\ Transform, IDFT)$，也记作 $\vec a=DFT_n^{-1}(\vec y)$，它是$DFT$的逆
• 这个问题实际上相当于是一个解线性方程组的问题，也就是给出了 $n$个线性方程
  $$\begin{equation*} \left\{ \begin{array}{ccccccccc} a_0(\omega_n^0)^{0}&+&\cdots&+&a_{n-2}(\omega_n^0)^{n-2}&+&+a_{n-1}(\omega_n^0)^{n-1}&=&A(\omega_n^0) \\ a_0(\omega_n^1)^{0}&+&\cdots&+&a_{n-2}(\omega_n^1)^{n-2}&+&+a_{n-1}(\omega_n^1)^{n-1}&=&A(\omega_n^1) \\ \vdots & & \vdots & &\vdots& & \vdots & & \vdots\\ a_0(\omega_n^{n-1})^{0}&+&\cdots&+&a_{n-2}(\omega_n^{n-1})^{n-2}&+&+a_{n-1}(\omega_n^{n-1})^{n-1}&=&A(\omega_n^{n-1}) \end{array} \right. \end{equation*}$$
  写成矩阵方程的形式就是
  $$\begin{equation} \begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & \cdots & (\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} \end{equation}$$
  记上面的系数矩阵为$V$，现在考虑下面这个矩阵$d_{ij}=ω^{−ij}_n$
  $$\begin{equation*} \mathbf D = \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \end{equation*}$$
  设它们相乘后的结果是$\mathbf E=\mathbf D \cdot \mathbf V$
  $$e_{ij} = \sum_{k=0}^{n-1} d_{ik} v_{kj} = \sum_{k=0}^{n-1} \omega_n^{-ik}\omega_n^{kj} =\sum_{k=0}^{n-1} \omega_n^{k(j-i)}$$
  当$i=j$时，$e_{ij}=n$
  当$i≠j$时，$e_{ij} = \displaystyle\sum_{k=0}^{n-1} (\omega_n^{j-i})^k = \frac{1-(\omega_n^{j-i})^n}{1-\omega_n^{j-i}} = 0$(根据等比数列求和公式)
  可得$\frac{1}{n}\mathbf E=\mathbf I_n$，所以$\mathbf E=\mathbf D \cdot \mathbf V\Rightarrow \frac{1}{n}\mathbf E=\frac{1}{n}\mathbf D \cdot \mathbf V\Rightarrow \mathbf I_n=\frac{1}{n}\mathbf D \cdot \mathbf V\Rightarrow \mathbf V^{-1}=\frac{1}{n}\mathbf D$
  将$\frac{1}{n}\mathbf D$即$\mathbf V^{-1}$在第一个矩阵两边左乘就会得到
  $$\begin{equation*} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \frac{1}{n} \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} \end{equation*}$$
• 据此，可以推导出$DFT^{−1}_N(\vec y)$，$a_k=\frac{1}{n}\displaystyle\sum_{i=0}^{n-1}{y_k\omega_n^{-ki}}$
  与$DFT_N(\vec a)$，$A(\omega_n^k) = \displaystyle\sum_{i=0}^{n-1}a_i\omega_n^{ki}$相比，我们发现两个式子十分相似
  只需要把$\vec a$和$\vec y$互换，用$ω_n^{−1}$替换$ω_n$，并将计算结果的每个元素除以$n$，可以在$O(nlogn)$的时间计算出$DFT^{−1}_n$
• 卷积定理 对任意两个长度为$n$的向量$\vec a$和$\vec b$，其中$n$是$2$的次幂，
  $$a\otimes b=DFT_{2n}^{-1}(DFT_{2n}(a)·DFT_{2n}(b))$$
  其中向量$\vec a$和$\vec b$用$0$填充，使其长度达到$2n$。并用$“·”$表示$2$个$2n$个元素组成的向量的点乘