FFT

给定两个多项式：

$$f(x)=1+2x+3x^2+4x^3,g(x)=4+3x+2x^2+x^3$$

如何求 $f(x)\cdot g(x)$？

方法一：暴力展开（朴素方法）

$$f(x)\cdot g(x)=(1+2x+3x^2+4x^3)(4+3x+2x^2+x^3)=4+11x+20x^2+30x^3+20x^4+11x^5+4x^6$$

• 时间复杂度：$O(n^2)$，其中 $n$ 是多项式的次数。
• 虽然简单直观，但当 $n$ 很大时（比如百万级），效率极低。

方法二：点值表示 + 插值（思路升级）

我们知道：一个 $d$ 次多项式可以由任意 $d+1$ 个不同点唯一确定。

如果我们能在若干点 $x_i$ 上快速计算出：$f(x_i)$ 和 $g(x_i)$，那么乘积多项式在这些点上的值就是：$h(x_i)=f(x_i)\cdot g(x_i)$

再通过插值得到 $h(x)$ 的系数形式即可！

但问题在于：求插值多项式的时间复杂度仍为 $O(n^2)$，能否优化？

分治思想 —— 奇偶次项分离

考虑将多项式 $f(x)$ 按奇偶次幂分离：

$$f(x)=1+2x+3x^2+4x^3=(1+3x^2)+x(2+4x^2)=f_e(x^2)+x\cdot f_o(x^2)$$

其中：

• $f_e(x)=1+3x$（偶次项）
• $f_o(x)=2+4x$（奇次项）

计算在 $x=1$ 和 $x=-1$ 处的值

$$f(1)=f_e(1^2)+1\cdot f_o(1^2)=f_e(1)+f_o(1)$$ $$f(-1)=f_e((-1{)}^2)+(-1)\cdot f_o((-1{)}^2)=f_e(1)-f_o(1)$$

只要计算一次 $f_e(1)$ 和 $f_o(1)$，就能同时得到 $f(1)$ 和 $f(-1)$ 两个点值

$f_e(x)=1+3x$ 和 $f_o(x)=2+4x$ 也可以用类似的思路继续划分

$$f_e(x)=1+3x=1+x\cdot 3$$ $$f_o(x)=2+4x=2+x\cdot 4$$

计算在 $x=1$ 和 $x=-1$ 处的值：

$$f_e(1)=1+(1)\cdot 3=4,f_e(-1)=1+(-1)\cdot 3=-2$$ $$f_o(1)=2+4\cdot 1=6,f_o(-1)=2+4\cdot (-1)=-2$$

我们知道 $f_e(1)$ 和 $f_o(1)$ 可以推出 $f(1)$ 和 $f(-1)$ 的值，那么 $f_e(-1)$ 和 $f_o(-1)$ 有什么用呢？

观察

$$f(x)=f_e(x^2)+x\cdot f_o(x^2)$$

可知，一个 $x^2$ 可以对应两个不同的 $x$ 值（比如 $x$ 和 $-x$），使得它们的平方相同：

那么有

$$f(a)=f_e(a^2)+a\cdot f_o(a^2)$$ $$f(-a)=f_e(a^2)-a\cdot f_o(a^2)$$

其中 $a^2=1$ 时 $a=\pm 1$，表示 $f_e(1)$ 和 $f_o(1)$ 可以推出 $f(1)$ 和 $f(-1)$ 的值

如果 $a^2=-1$ 呢?

如果你学过复数，会立马想到：

$$i^2=-1$$

也就有

$$f(i)=f_e(i^2)+i\cdot f_o(i^2)=f_e(-1)+i\cdot f_o(-1)=-2+i\cdot (-2)=-2-2i$$ $$f(-i)=f_e((-i{)}^2)+(-i)\cdot f_o((-i{)}^2)=f_e(-1)+(-i)\cdot f_o(-1)=-2+(-i)(-2)=-2+2i$$

利用虚数单位 $i$ 进行取值，由此我们可以得到四个点：

$$f(1)=10$$ $$f(i)=-2-2i$$ $$f(-1)=-2$$ $$f(-i)=-2+2i$$

下面使用一个动画展示上述计算过程

（此处插入动画演示）

在动画中可以看到，我们找到了平方为 $1$ 的两个数 $1$ 和 $-1$ 和 平方为 $-1$ 的两个数 $i$ 和 $-i$，这看起来是个递归的过程

（插入动画单位根树状图）

根据欧拉公式 $e^{i\theta }=\cos \theta +i\sin \theta$，可以把这些值放到单位圆上，表示为 $e^{\frac{k}{n}2\pi i}$

这些点被称为 $n$ 次单位根

观察到 ${\omega }_8^2=e^{\frac{2}{8}2\pi i}=e^{\frac{1}{4}2\pi i}={\omega }_4^1{\omega }_8^{1+4}=-{\omega }_8^1$

令 ${\omega }_n=e^{\frac{1}{n}2\pi i}$，单位圆每个数分别为 ${\omega }_n^k$

这些点还有一些性质

$${\omega }_{mn}^{mk}={\omega }_n^k,{\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$$

结合这些性质我们可以带入到公式 $f(x)=f_e(x^2)+x\cdot f_o(x^2)$ 中

$$f({\omega }_n^k)=f_e({\omega }_n^{2k})+{\omega }_n^k\cdot f_o({\omega }_n^{2k})=f_e({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\cdot f_o({\omega }_{\frac{n}{2}}^k)$$ $$f({\omega }_n^{k+\frac{n}{2}})=f_e(({\omega }_n^{k+\frac{n}{2}}{)}^2)+{\omega }_n^{k+\frac{n}{2}}\cdot f_o(({\omega }_n^{k+\frac{n}{2}}{)}^2)=f_e({\omega }_n^{2k})+{\omega }_n^{k+\frac{n}{2}}\cdot f_o({\omega }_n^{2k})=f_e({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\cdot f_o({\omega }_{\frac{n}{2}}^k)$$

至此，可以使用 python 实现 FFT 计算过程

def FFT(P):
    n = len(P)
    if n == 1:
        return P
    Pe = P[::2]
    Po = P[1::2]
    Fe = FFT(Pe)
    Fo = FFT(Po)
    wn = np.exp(2 * np.pi / n * 1j)
    F = [0] * n
    for k in range(n // 2):
        F[k] = Fe[k] + wn**k * Fo[k]
        F[k + n // 2] = Fe[k] - wn**k * Fo[k]
    return F

回到正题，我们已经得到了 $n$ 个插值点，现在如何还原为多项式呢？

$$P(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$$

将这些单位根 ${\omega }_n^k$ 带入 $P(x)$

$$\left[\begin{array}{c}P(1)\\ P({\omega }_n)\\ P({\omega }_n^2)\\ ⋮\\ P({\omega }_n^{n-1})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right]$$

也就是

$$P(w_n^k)=\sum _{j=0}^{n-1}{a}_j{w}_n^{jk}=\sum _{j=0}^{n-1}{a}_j{e}^{\frac{jk}{n}2\pi i}$$

这是傅里叶变换

$$F(\omega )=\mathbb{F}[f(t)]={\int }_{-\infty }^{\infty }f(t){\mathrm{e}}^{-\mathrm{i}\omega t}dt$$

的离散形式，称为离散傅里叶变换（Discrete Fourier transform，DFT），这也是本算法称为快速傅立叶变换（Fast Fourier Transform，FFT）原因

其中

$$\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]$$

是范德蒙德矩阵，取逆矩阵可得

$$\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right]=\frac{1}{n}\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n^{-1}& {\omega }_n^{-2}& \cdots & {\omega }_n^{-(n-1)}\\ 1& {\omega }_n^{-2}& {\omega }_n^{-4}& \cdots & {\omega }_n^{-2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{-(n-1)}& {\omega }_n^{-2(n-1)}& \cdots & {\omega }_n^{-(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}P(1)\\ P({\omega }_n)\\ P({\omega }_n^2)\\ ⋮\\ P({\omega }_n^{n-1})\end{array}\right]$$

只需将单位元 ${\omega }_n$ 取逆（倒数） ${\omega }_n^{-1}$ ，即可得到 FFT 逆变换

def FFT(P, inv=1): # inv=-1表示逆变换
    n = len(P)
    if len(P) == 1:
        return P
    Pe = P[::2]
    Po = P[1::2]
    Fe = FFT(Pe, inv)
    Fo = FFT(Po, inv)
    wn = np.exp(inv * 2 * np.pi / n * 1j)
    F = [0] * n
    for k in range(n // 2):
        F[k] = Fe[k] + wn**k * Fo[k]
        F[k + n // 2] = Fe[k] - wn**k * Fo[k]
    return F
    
def IFFT(P: list):
    n = len(P)
    F = FFT(P, -1)
    return [x / n for x in F]

回到最开始的问题，如何求解两个多项式的乘积

def multiply(A, B):
    n = len(A) + len(B) - 1
    nn = 2 << int(math.log2(n - 1))
    f = FFT(A + [0j] * (nn - len(A)))
    g = FFT(B + [0j] * (nn - len(B)))
    h = [x * y for x, y in zip(f, g)]
    S = IFFT(h)
    return [int(x.real + 0.49) for x in S[:n]]