-1. 前置知识
基础的复数知识。
0. 什么是多项式乘法
众所周知,多项式本质是一种特殊的函数,可以表示为自变量的若干次幂之和,即
\[F(x)=\sum_{i=0}c_i\cdot x^i \]其中 \(c_i\) 被称为 \(x^i\) 的系数。
已知 \(F,G\) 是两个多项式函数,考虑定义一个新的函数 \(H(x)=F(x)G(x)\)。我们称 \(H\) 为 \(F\) 与 \(G\) 的乘积,记作 \(H=F\cdot G\),这种通过已知的两个多项式函数,以乘积形式生成另一个函数的运算,称为多项式乘法。
本文讨论的就是计算多项式乘法,即求出 \(H\) 解析式的过程。
根据定义,有
\[\begin{aligned} H(x)&=F(x)\cdot G(x) \\ &= (\sum_{i=0} a_ix^i)(\sum_{j=0}b_jx^j) \\ &= \sum_{i=0}(a_ix^i\cdot \sum_{j=0}b_jx^j)\\ &= \sum_{i=0}\sum_{j=0} a_ib_jx^{i+j}\\ \end{aligned} \]我们发现,两个多项式的乘积仍是多项式。
由于多项式的特性,我们只需知道该多项式的每一项系数即可,即使用一个\(n\) 维向量去描述一个\(n-1\)次多项式。
显然,有一种暴力的方法 \(H[i]=\sum_{j=0}^iF[j]* G[i-j]\)。我们称这种形如\(c_n=\sum_{i\otimes j=n}a_i\cdot b_j\)运算为卷积。(其实只是卷积的一种,即加法卷积)
显然,这样做复杂度是 \(\mathcal{O}(n^2)\)的,实在太慢了。
\(\frac{1}{2}\). 多项式的点值表示法
既然多项式是一个函数,那我么可以画出它的函数图象,并在上面取几个点。
显然,由待定系数法可知,给定 \(n+1\) 个不同的点,可以唯一确定一个次数不超过 \(n\) 的多项式。
我们可以考虑用点值反推解析式,即如果已知函数 \(F,G\) 的\(n\) 个对应点,由定义式 \(H(x)=F(x)G(x)\) 可以直接 \(\mathcal O(n)\) 算出 \(H\) 对应的点!
但如果随意带入点,用待定系数法高斯消元暴力反推解析式,复杂度 \(\mathcal O(n^3)\),直接爆炸。
当然我们可以用拉格朗日插值公式做到 \(\mathcal O(n^2)\),然并卵,还不如暴力乘法……
\(\frac{9}{10}\). 单位根的引入及其性质
既然带一些普通的东西进去没啥用,那我们直接躺平带一些奇怪的东西进去,比如复数。考虑带入 \(n\) 次单位根。
定义:方程 \(x^n-1=0\) 在复数域的全部解称为 \(n\) 次单位根。由代数基本定理知,\(n\) 次单位根有 \(n\) 个,分别记为 \(\omega_n^0 \cdots \omega_n^{n-1}\) 。由复数乘法模相乘,幅角相加可知,\(n\) 次单位根全部在单位圆上,且幅角为 \(\frac{2\pi}{n}\) 的整数倍。显然,\(\omega_n^k=\cos \frac{2\pi k}{n}+i\sin \frac{2\pi k}{n}\)。当然,其它单位根是 \(\omega_n^1\) 的整数次幂,因此我们称其为主\(n\)次单位根,简记为 \(\omega_n\)。
根据欧拉公式 \(e^{ix}=\cos x+i\sin x\),有\(\omega_n^k=e^{\frac{2\pi k}{n}i}\)
一些性质:
-
\(\omega_n^k=(\omega_n^1)^k\)
-
\(\omega_n^i\cdot\omega_n^j=\omega_n^{i+j}\)
-
\(\omega_{dn}^{dk}=\omega_n^k\) (消去引理)
证明:\(\large\omega_{dn}^{dk}=e^{\frac{2\pi dk}{dn}i}=e^{\frac{2\pi k}{n}i}=\omega_n^k\)
- 若 \(2|n\),\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)
证明:\(\omega_n^{k+\frac{n}{2}}=\omega_n^{k}\omega_n^{\frac{n}{2}}=\omega^1_2\omega_n^k=-\omega_n^k\)
推论:\((\omega_n^k)^2=(\omega_n^{k+\frac{n}{2}})^2=\omega_{n/2}^k\)(折半引理)
- 对于任意不为 \(n\) 倍数的 \(k\),有:
证明:
\[\sum_{i=0}^{n-1}(\omega_n^k)^i=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{1^k-1}{\omega_{n}^k-1}=0 \]知道这些结论,就可以学习 FFT 了!
1. DFT
考虑求出一个多项式在 \(\omega_n^0\cdots\omega_n^{n-1}\) 的值,我们称这种运算为 DFT。
首先,我们将该多项式的次数扩充为 \(2\) 的幂次。
将 \(A(x)\) 写成系数向量形式 \(a=[a_0,a_1,a_2,\cdots,a_{n-1}]\)
考虑将其按奇数项和偶数项分为
\[a^{[0]}=[a_0,a_2,a_4,\cdots,a_{n-2}] \]\[a^{[1]}=[a_1,a_3,a_5,\cdots,a_{n-1}] \]分别记其对应的多项式为 \(A^{[0]}(x)\) 和 \(A^{[1]}(x)\)。
\[A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} \]\[A^{[0]}(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{n/2-1} \]\[A^{[1]}(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{n/2-1} \]换成 \(x^2\):
\[A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} \]\[A^{[0]}(x^2)=a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2} \]\[A^{[1]}(x^2)=a_1+a_3x^2+a_5x^4+\cdots+a_{n-1}x^{n-2} \]\(A^{[1]}\) 乘上 \(x\):
\[A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} \]\[A^{[0]}(x^2)=a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2} \]\[xA^{[1]}(x^2)=a_1x+a_3x^3+a_5x^5+\cdots+a_{n-1}x^{n-1} \]由此得到 \(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)
说人话:如果知道 \(A^{[0]}(x^2)\) 和 \(A^{[1]}(x^2)\) ,就可以求出 \(A(x)\)。
问题转化为求出所有 \(A^{[0]}((\omega_n^k)^2)\) 和 \(A^{[1]}((\omega_n^k)^2)\)。
由折半引理知,求出所有 \(A^{[0]}((\omega_n^k)^2)\) 和 \(A^{[1]}((\omega_n^k)^2)\) 就是求出所有 \(A^{[0]}(\omega_{n/2}^k)\) 和 \(A^{[1]}(\omega_{n/2}^k)\)。
注意到这是原问题的子问题,可以直接分治处理。递归边界为 \(n=1\),此时对应的多项式为常量,直接返回即可。
考虑合并。
\(A(\omega_n^k)=A^{[0]}(\omega_{n/2}^k)+\omega_n^kA^{[1]}(\omega_{n/2}^k)\)
当然,如果 \(k\geq \frac{n}{2}\),我们让 \(k\bmod \frac{n}{2}\) 即可。
显然,合并复杂度为 \(\mathcal O(n)\)。
则 \(FFT\) 的复杂度 \(T(n)=2T(\frac{n}{2})+\mathcal O(n)\),故 \(T(n)=\mathcal O(n\log n)\)
全剧终。
我们发现,我们得到的只是一些点值罢了,根本不是系数表示。
2. IDFT
即 DFT 的逆运算。
如何将点值转换为系数?注意到上文提到 DFT 时,说的是将 \(A(x)\) 看成系数向量,即一维矩阵。
如果我们将得到的点值记为一维向量 \(b\),则有:\(a=IDFT(b)\),即 \(b=DFT(a)\)。
我们现在知道 \(b\) 求 \(a\)。
由定义知 \(b[k]=\sum_{i=0}^{n-1}(\omega_n^k)^i\cdot a[i]\)
如果你会单位根反演,你可能能够推出 \(a\) 关于 \(b\) 的式子。这里给出单位根反演的原理。
还记得求和引理吗?
\[\sum_{i=0}^{n-1}(\omega_n^k)^i=0 \]对于任意不为 \(n\) 倍数的 \(k\),有:
此外,若 \(k\) 为 \(n\) 的倍数,显然上个式子的值为 \(n\)。
即
\[\frac{1}{n}\sum_{i=0}^{n-1}(\omega_n^k)^i=[n|k] \]考虑如下式子:
\[a_k=\sum_{i=0}^{n-1} a_i\cdot[k=i] \]感觉是句废话……
考虑把 \([k=i]\) 往 \([n|k]\) 上靠。
我们发现,\(i,k\) 的取值范围是 \([0,n-1]\),即 \(k=k\bmod n,i=i\bmod n\)
那么,\([k=i]\) 可以写成 $[k\bmod n=i\bmod n] $,即 $ [k\equiv i\pmod n]$ 或 \([n|(k-i)]\)
带回原式试试!
\[\begin{aligned} a_k &=\sum_{i=0}^{n-1}a_i\cdot[n|(k-i)]\\ &=\sum_{i=0}^{n-1}a_i\cdot\frac{1}{n}\sum_{j=0}^{n-1}(\omega_n^{k-i})^j \\ &=\frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} \omega_n^{ik}\omega_n^{-ij}a_i \\ &=\frac{1}{n}\sum_{i=0}^{n-1}(\omega_n^{ki})\cdot \sum_{j=0}^{n-1}\omega_n^{-ij}a_i \end{aligned} \]完蛋了……推不出来了
但 \([n|(k-i)]\) 与 \([n|(i-k)]\) 是等价的。
再试一次。
\[\begin{aligned} a_k &=\sum_{i=0}^{n-1}a_i\cdot[n|(i-k)]\\ &=\sum_{i=0}^{n-1}a_i\cdot\frac{1}{n}\sum_{j=0}^{n-1}(\omega_n^{i-k})^j \\ &=\frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} \omega_n^{-ik}\omega_n^{ij}a_i \\ &=\frac{1}{n}\sum_{i=0}^{n-1}(\omega_n^{-ki})\cdot \sum_{j=0}^{n-1}\omega_n^{ij}a_i \\ &=\frac{1}{n}\sum_{i=0}^{n-1}(\omega_n^{-k})^i\cdot b_i \end{aligned} \]这次成功了。变形一下:
\[n\cdot a[k]=\sum_{i=0}^{n-1}(\omega_n^{-k})^ib[i] \]我们发现右边和 DFT 过程几乎完全一致,只是将 \(\omega_n^k\) 改为 \(\omega_n^{-k}\)。
到这里,FFT 的理论部分就讲完了。