多项式、快速傅里叶变换、FFT

皓码
发布于 2天前
阅读 77

快速傅里叶变换（FFT）快速傅里叶变换（FastFourierTransform，简称FFT）的高效算法主要由詹姆斯・库利（JamesW.Cooley）和约翰・图基（JohnW.Tukey）在1965年提出。

## 快速傅里叶变换（FFT）
快速傅里叶变换（Fast Fourier Transform，简称 FFT）的高效算法主要由詹姆斯・库利（James W. Cooley）和约翰・图基（John W. Tukey）在 1965 年提出。他们在论文《机器计算傅里叶级数的一种算法》中，首次系统地阐述了通过分治策略将离散傅里叶变换（DFT）的计算复杂度从 $O(N^2)$ 降低至 $O(N\log N)$ 的方法，这一突破极大推动了傅里叶变换在信号处理、图像处理、密码学等领域的实际应用。

多项式与快速傅里叶变换（FFT）的核心关联在于：FFT 是高效计算多项式乘法（卷积）的关键工具。通过 FFT，多项式乘法的时间复杂度可从传统的 $O(n^2)$ 优化至 $O(n \log n)$，这对数值计算、信号处理、组合数学等领域具有革命性意义。

话不多说直接上干货

### 欧拉公式
$e^{i\theta}=cos\theta+i*sin\theta$，$i$是虚数单位

### 参数设置与定义
考虑下面两个场景：
1. 已知两个多项式快速的求出在 N 个不同点的乘积，即 FFT 过程。
2. 已知多项式在 N 个不同点的乘积，求这个多项式的系数，是FFT的逆操作，即 IFFT 过程。使用拉格朗日插值法去做效率不高，因为求范德蒙德矩阵的逆的复杂度至少为 $O(n^2)$ 。

假设有两个多项式 $A(x),B(x)$ , 快速求两个多项式的乘积（系数表示）。首先做一个预处理，将两个多项式的系数个数进行补全，补到系数个数为2的幂次这种形式。不妨假设这两个多项式的个数分别为 $n_1,n_2$ , 那么这两个多项式的乘积的系数个数为 $n=(n_1-1)+(n_2-1)+1=n_1+n_2-1$ 。取一个离$n$最近的且是 2 的幂次的整数，设为 $N=2^t$ , 然后将 $A(x),B(x)$ 的高次项系数补 0 。

离散傅里叶变换 (DFT) 是来计算多项式在 $N$ 个特殊点(单位根)的值。而快速傅里叶变换 (FFT) 是一种快速有效率的对 DFT 的实现。FFT 加速多项式乘法，其基本思想是将两个多项式的系数表示通过 FFT 转化为特殊点处的点值表示，然后计算两个多项式点值表示的乘积得到原多项式卷积的点值表示，再将多项式卷积的点值表示进行逆离散傅里叶变换 (IDFT) 就得到了乘积多项式的系数表示。

看不明白的可以跟着公式，一起推导一遍

### 一般方法
假设已经求出两个多项式在 $N$ 个不同点的值 , 设 $C(x)=A(x)*B(x)$ ，进而计算出 $C(x)$ 在 $N$ 个点的值。
$$
\begin{aligned}
    A(x_i)=a_i,i=0,1,...,N-1\\
    B(x_i)=b_i,i=0,1,...,N-1
\end{aligned}
$$
计算出 $c_i=C(x_i)=a_i*b_i,i=0,1,...,N-1$ 。然后根据范德蒙德矩阵求出 $C(x)$ 的系数 $C_0,...,C_{N-1}$ 。
$$
D=\begin{pmatrix}
1 & x_0&x_0^2&\dots&x_0^{N-1}\\
1 & x_1&x_1^2&\dots&x_1^{N-1}\\
\dots & \dots \\
1 & x_{N-1}&x_{N-1}^2&\dots&x_{N-1}^{N-1}\\
\end{pmatrix}\\

$$
$$
\begin{aligned}
    
\vec{C}&=[C_0,C_1,\dots,C_{N-1}]^T\\
\vec{c}&=[c_0,c_1,\dots,c_{N-1}]^T\\
\vec{c}&=D*\vec{C}\\

\vec{C}&=D^{-1}*\vec{c}
\end{aligned}
$$

上面的推导过程就是两个多项式相乘得到的新多项式的系数的求解过程，求解范德蒙德矩阵的逆矩阵效率是不高的，只做理论分析。

### 快速方法
这个快速方法就是 FFT，在复数域内考虑方程 $x^N=1$ ，由代数学知识我们知道，它有 $N$ 个根，分别记为

$\omega_N^{k}=e^{i\frac{2k\pi}{N}},k=0,\dots,N-1$，当然欧拉公式可以展开。可以考虑计算多项式在这些特殊的点的值。
验证
$$
\begin{aligned}
    (\omega_N^k)^N&=(e^{i\frac{2k\pi}{N}})^N\\
    &=e^{i*2k\pi}\\
    &=cos2k\pi + i*sin2k\pi\\
    &=1
\end{aligned}
$$
所以上述形式 $\omega_N^{k}$ 是方程 $x^N=1$ 的 $N$ 个不同的根，如果已知多项式，那么就可以计算在这 $N$ 个点处的值。

重要性质：
$$
\begin{aligned}
    \omega_N^k&=\omega_N^{k-1}*\omega_N^{1} \tag{1}
\end{aligned}
$$

$$
\begin{aligned}
    \omega_N^0&=\omega_N^{N}=1 \tag{2}
\end{aligned}
$$

$$
\begin{aligned}                 \omega_{2N}^{2k}&=\omega_N^k=\omega_{mN}^{mk} \tag{3}
\end{aligned}
$$

$$
\begin{aligned}
    \omega_N^{k+\frac{N}{2}}&=-\omega_N^k ,k=0,...,\frac{N}{2}-1 \tag{4}
\end{aligned}
$$

$$
\begin{aligned}
    \omega_N^{2k}&=-\omega_{\frac{N}{2}}^k \tag{5}
\end{aligned}
$$

$$
\begin{aligned}       
    \sum_{j=0}^{N-1}{\omega_N^{jk}}&=\left\{
        \begin{aligned}
            0,k\text{ mod }N\neq0\\
            n,k\text{ mod }N=0\\            
        \end{aligned}
    \right.\text{等比数列可证}\tag{6}
\end{aligned}
$$
根据欧拉公式，上述几个性质非常容易推导出来，这里不再赘述。

假设如下多项式的次数 $N$ 是 2 的幂次，即存在 $t$，满足 $N=2^t$。

$$
\begin{aligned}
f(x)&=a_0+a_1x+a_2x^2+...+a_{N-1}x^{N-1}\\
&=(a_0+a_2x^2+\dots+a_{N-2}x^{N-2})+x(a_1+a_3x^2+\dots+a_{N-1}x^{N-2})\\
&=f_0(x^2)+xf_1(x^2) \tag{7}
\end{aligned}
$$
上述公式给出了一种迭代形式，将单位根代入到公式 $(7)$，迭代公式如下：
$$
\begin{aligned}
f(\omega_N^k)&=f_0(\omega_N^{2k})+\omega_N^kf_1(\omega_N^{2k})\\
&=f_0(\omega_{\frac{N}{2}}^k)+\omega_N^kf_1(\omega_{\frac{N}{2}}^k) 
\tag{8}
\end{aligned}
$$

$$
\begin{aligned}    
f(\omega_N^{k+\frac{N}{2}})&=f(-\omega_N^k)\\
&=f_0(\omega_{\frac{N}{2}}^k)-\omega_N^kf_1(\omega_{\frac{N}{2}}^k)\\
k&=0,1,\dots,\frac{N}{2}-1 \tag{9}
\end{aligned}
$$
继续进行下去，根据公式 $(8)$，有
$$
\begin{aligned}
    f_0(\omega_{\frac{N}{2}}^k)=f_{00}(\omega_{\frac{N}{4}}^k)+\omega_{\frac{N}{2}}^kf_{01}(\omega_{\frac{N}{4}}^k)
\end{aligned}
$$
$$
\begin{aligned}
    f_1(\omega_{\frac{N}{2}}^k)=f_{10}(\omega_{\frac{N}{4}}^k)-\omega_{\frac{N}{2}}^kf_{11}(\omega_{\frac{N}{4}}^k)
\end{aligned}
$$

这种形式很像一颗二叉树，如果不理解自己可以假设 $N==8$，按照上述方法逐渐写出所有的裂解形式，就很容易理解了。
### 快速多项式乘积 (IFFT)
有了上述算法，就很容易快速的计算出两个多项式 $A(x),B(x)$ 分别在 $N$ 个单位根处的值。
$$
A(\omega_N^k)=a_k,k=0,1,\dots,N-1
$$
$$
B(\omega_N^k)=b_k,k=0,1,\dots,N-1
$$

所以有，
$$
C(\omega_N^k)=a_k*b_k=c_k,k=0,1,\dots,N-1
$$

新的范德蒙德行列式如下，

$$
W=\begin{pmatrix}
    1 & 1&1 & \dots & 1\\
    1 & \omega_N^{1*1}& \omega_N^{1*2}&\dots & \omega_N^{1*(N-1)}\\
    1 & \omega_N^{2*1}&\omega_N^{2*2}&\dots&\omega_N^{2*(N-1)}\\
    \dots & \dots \\
    1 & \omega_N^{(N-1)*1}&\omega_N^{(N-1)*2}& \dots &\omega_N^{(N-1)*(N-1)}\\
\end{pmatrix}
$$
对于一般的范德蒙德矩阵的逆的复杂度是$O(n^3)$
对于这些特殊点构成的范德蒙的矩阵逆是非常容易求出了。逆的形式如下：
$$
W^{-1}=\frac{1}{N}\begin{pmatrix}
    1 & 1&1 & \dots & 1\\
    1 & \omega_N^{-1*1}& \omega_N^{-1*2}&\dots & \omega_N^{-1*(N-1)}\\
    1 & \omega_N^{-2*1}&\omega_N^{-2*2}&\dots&\omega_N^{-2*(N-1)}\\
    \dots & \dots \\
    1 & \omega_N^{-(N-1)*1}&\omega_N^{-(N-1)*2}& \dots &\omega_N^{-(N-1)*(N-1)}\\
\end{pmatrix}
$$

于是多项式$C(x)$的系数向量 $\vec{C}$ 就是
$$
\begin{aligned}
    \vec{C}&=W^{-1}*\vec{c}\\
    \vec{c}&=[c_0,\dots,c_{N-1}]^T
\end{aligned}
$$

到这里是不是有点感觉了，是不是能和线性代数里讲的正交基的概念有联系？自己可以思考下。
### NTT 
有限域内的快速多项式乘法计算。在有限域里我们可以发现原根具有着一些和单位根很相似的性质。原根如下：
$g$是模$p$(素数)的原根当且仅当$g$的阶是$\phi(p)=p-1$,事实上，$g$是乘法循环群$Z_p^*$的生成元。
$$
g_N^k=g^{\frac{p-1}{N}k}，k=0,1,\dots,N-1
$$

很容易验证单位根有的性质，有限域中原根同样也有。所以NTT是有限域下的快速计算多项式快速的算法。单位根替换成原根。

1. 原根与模数条件
- 原根：设 $p$ 为质数，若存在整数 $g$ 使得 $g$ 的阶为 $p-1$（即 $g^{p-1} \equiv 1 \mod p$，且对任意 $0 < k < p-1$，$g^k \not\equiv 1 \mod p$），则称 g 是模 p 的原根。

- 模数要求：通常选择费马素数 $p = 2^k \cdot m + 1$（如 $p = 2^{23} \cdot 119 + 1 = 998244353$），使得 $n = 2^k$ 是 $p-1$ 的因数，此时原根 g 的 n 次单位根为 $g^{(p-1)/n} \mod p$。

至此快速傅里叶变换已经写完。

快速傅里叶变换（FFT）

快速傅里叶变换（Fast Fourier Transform，简称 FFT）的高效算法主要由詹姆斯・库利（James W. Cooley）和约翰・图基（John W. Tukey）在 1965 年提出。他们在论文《机器计算傅里叶级数的一种算法》中，首次系统地阐述了通过分治策略将离散傅里叶变换（DFT）的计算复杂度从 $O(N^2)$ 降低至 $O(N\log N)$ 的方法，这一突破极大推动了傅里叶变换在信号处理、图像处理、密码学等领域的实际应用。

多项式与快速傅里叶变换（FFT）的核心关联在于：FFT 是高效计算多项式乘法（卷积）的关键工具。通过 FFT，多项式乘法的时间复杂度可从传统的 $O(n^2)$ 优化至 $O(n \log n)$，这对数值计算、信号处理、组合数学等领域具有革命性意义。

话不多说直接上干货

欧拉公式

$e^{i\theta}=cos\theta+i*sin\theta$，$i$是虚数单位

参数设置与定义

考虑下面两个场景：

已知两个多项式快速的求出在 N 个不同点的乘积，即 FFT 过程。
已知多项式在 N 个不同点的乘积，求这个多项式的系数，是FFT的逆操作，即 IFFT 过程。使用拉格朗日插值法去做效率不高，因为求范德蒙德矩阵的逆的复杂度至少为 $O(n^2)$ 。

看不明白的可以跟着公式，一起推导一遍

一般方法

假设已经求出两个多项式在 $N$ 个不同点的值 , 设 $C(x)=A(x)B(x)$ ，进而计算出 $C(x)$ 在 $N$ 个点的值。 $$ \begin{aligned} A(x_i)=a_i,i=0,1,...,N-1\ B(x_i)=b_i,i=0,1,...,N-1 \end{aligned} $$ 计算出 $c_i=C(x_i)=a_ib_i,i=0,1,...,N-1$ 。然后根据范德蒙德矩阵求出 $C(x)$ 的系数 $C0,...,C{N-1}$ 。 $$ D=\begin{pmatrix} 1 & x_0&x_0^2&\dots&x_0^{N-1}\ 1 & x_1&x_1^2&\dots&x1^{N-1}\ \dots & \dots \ 1 & x{N-1}&x{N-1}^2&\dots&x{N-1}^{N-1}\ \end{pmatrix}\

$$ $$ \begin{aligned}

\vec{C}&=[C_0,C1,\dots,C{N-1}]^T\ \vec{c}&=[c_0,c1,\dots,c{N-1}]^T\ \vec{c}&=D*\vec{C}\

\vec{C}&=D^{-1}*\vec{c} \end{aligned} $$

上面的推导过程就是两个多项式相乘得到的新多项式的系数的求解过程，求解范德蒙德矩阵的逆矩阵效率是不高的，只做理论分析。

快速方法

这个快速方法就是 FFT，在复数域内考虑方程 $x^N=1$ ，由代数学知识我们知道，它有 $N$ 个根，分别记为

$\omega_N^{k}=e^{i\frac{2k\pi}{N}},k=0,\dots,N-1$，当然欧拉公式可以展开。可以考虑计算多项式在这些特殊的点的值。验证 $$ \begin{aligned} (\omega_N^k)^N&=(e^{i\frac{2k\pi}{N}})^N\ &=e^{i2k\pi}\ &=cos2k\pi + isin2k\pi\ &=1 \end{aligned} $$ 所以上述形式 $\omega_N^{k}$ 是方程 $x^N=1$ 的 $N$ 个不同的根，如果已知多项式，那么就可以计算在这 $N$ 个点处的值。

重要性质： $$ \begin{aligned} \omega_N^k&=\omega_N^{k-1}*\omega_N^{1} \tag{1} \end{aligned} $$

$$ \begin{aligned} \omega_N^0&=\omega_N^{N}=1 \tag{2} \end{aligned} $$

$$ \begin{aligned} \omega_{2N}^{2k}&=\omegaN^k=\omega{mN}^{mk} \tag{3} \end{aligned} $$

$$ \begin{aligned} \omega_N^{k+\frac{N}{2}}&=-\omega_N^k ,k=0,...,\frac{N}{2}-1 \tag{4} \end{aligned} $$

$$ \begin{aligned} \omegaN^{2k}&=-\omega{\frac{N}{2}}^k \tag{5} \end{aligned} $$

$$ \begin{aligned}
\sum_{j=0}^{N-1}{\omega_N^{jk}}&=\left{ \begin{aligned} 0,k\text{ mod }N\neq0\ n,k\text{ mod }N=0\
\end{aligned} \right.\text{等比数列可证}\tag{6} \end{aligned} $$ 根据欧拉公式，上述几个性质非常容易推导出来，这里不再赘述。

假设如下多项式的次数 $N$ 是 2 的幂次，即存在 $t$，满足 $N=2^t$。

$$ \begin{aligned} f(x)&=a_0+a_1x+a2x^2+...+a{N-1}x^{N-1}\ &=(a_0+a2x^2+\dots+a{N-2}x^{N-2})+x(a_1+a3x^2+\dots+a{N-1}x^{N-2})\ &=f_0(x^2)+xf_1(x^2) \tag{7} \end{aligned} $$ 上述公式给出了一种迭代形式，将单位根代入到公式 $(7)$，迭代公式如下： $$ \begin{aligned} f(\omega_N^k)&=f_0(\omega_N^{2k})+\omega_N^kf_1(\omega_N^{2k})\ &=f0(\omega{\frac{N}{2}}^k)+\omega_N^kf1(\omega{\frac{N}{2}}^k) \tag{8} \end{aligned} $$

$$ \begin{aligned}
f(\omega_N^{k+\frac{N}{2}})&=f(-\omega_N^k)\ &=f0(\omega{\frac{N}{2}}^k)-\omega_N^kf1(\omega{\frac{N}{2}}^k)\ k&=0,1,\dots,\frac{N}{2}-1 \tag{9} \end{aligned} $$ 继续进行下去，根据公式 $(8)$，有 $$ \begin{aligned} f0(\omega{\frac{N}{2}}^k)=f{00}(\omega{\frac{N}{4}}^k)+\omega{\frac{N}{2}}^kf{01}(\omega_{\frac{N}{4}}^k) \end{aligned} $$ $$ \begin{aligned} f1(\omega{\frac{N}{2}}^k)=f{10}(\omega{\frac{N}{4}}^k)-\omega{\frac{N}{2}}^kf{11}(\omega_{\frac{N}{4}}^k) \end{aligned} $$

这种形式很像一颗二叉树，如果不理解自己可以假设 $N==8$，按照上述方法逐渐写出所有的裂解形式，就很容易理解了。

快速多项式乘积 (IFFT)

有了上述算法，就很容易快速的计算出两个多项式 $A(x),B(x)$ 分别在 $N$ 个单位根处的值。 $$ A(\omega_N^k)=a_k,k=0,1,\dots,N-1 $$ $$ B(\omega_N^k)=b_k,k=0,1,\dots,N-1 $$

所以有， $$ C(\omega_N^k)=a_k*b_k=c_k,k=0,1,\dots,N-1 $$

新的范德蒙德行列式如下，

$$ W=\begin{pmatrix} 1 & 1&1 & \dots & 1\ 1 & \omega_N^{11}& \omega_N^{12}&\dots & \omega_N^{1(N-1)}\ 1 & \omega_N^{21}&\omega_N^{22}&\dots&\omega_N^{2(N-1)}\ \dots & \dots \ 1 & \omega_N^{(N-1)1}&\omega_N^{(N-1)2}& \dots &\omega_N^{(N-1)(N-1)}\ \end{pmatrix} $$ 对于一般的范德蒙德矩阵的逆的复杂度是$O(n^3)$ 对于这些特殊点构成的范德蒙的矩阵逆是非常容易求出了。逆的形式如下： $$ W^{-1}=\frac{1}{N}\begin{pmatrix} 1 & 1&1 & \dots & 1\ 1 & \omega_N^{-11}& \omega_N^{-12}&\dots & \omega_N^{-1(N-1)}\ 1 & \omega_N^{-21}&\omega_N^{-22}&\dots&\omega_N^{-2(N-1)}\ \dots & \dots \ 1 & \omega_N^{-(N-1)1}&\omega_N^{-(N-1)2}& \dots &\omega_N^{-(N-1)(N-1)}\ \end{pmatrix} $$

于是多项式$C(x)$的系数向量 $\vec{C}$ 就是 $$ \begin{aligned} \vec{C}&=W^{-1}*\vec{c}\ \vec{c}&=[c0,\dots,c{N-1}]^T \end{aligned} $$

到这里是不是有点感觉了，是不是能和线性代数里讲的正交基的概念有联系？自己可以思考下。

NTT

有限域内的快速多项式乘法计算。在有限域里我们可以发现原根具有着一些和单位根很相似的性质。原根如下： $g$是模$p$(素数)的原根当且仅当$g$的阶是$\phi(p)=p-1$,事实上，$g$是乘法循环群$Z_p^*$的生成元。 $$ g_N^k=g^{\frac{p-1}{N}k}，k=0,1,\dots,N-1 $$