多项式相乘本质上是计算系数的卷积，但同时也是值域上的点积，因此可以通过离散傅里叶变换将卷积转化为点积， 同时由于该傅里叶变换和逆变换都可以通过递归分治求解，每次分治都能将问题规模缩小一半，使复杂度从 $\mathcal{O}(N^2)$ 降到了 $\mathcal{O}(N\log N)$.

1. 基础

首先, 对于一个 $n$ 阶多项式 $f\left( x \right) =a_0+a_1x+a_2x^2+\cdots +a_nx^n$, 我们只需取 $n+1$ 个函数值, 就能反推出 $f(x)$ 的系数。

首先我们将其写成矩阵形式:

\[\left( \begin{matrix} 1& x_0& {x_0}^2& \cdots& {x_0}^n\\ 1& x_1& {x_1}^2& \cdots& {x_1}^n\\ 1& x_2& {x_2}^2& \cdots& {x_2}^n\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_n& {x_n}^2& \cdots& {x_n}^n\\ \end{matrix} \right) \left( \begin{array}{c} a_0\\ a_1\\ a_2\\ \vdots\\ a_n\\ \end{array} \right) =\left( \begin{array}{c} y_0\\ y_1\\ y_2\\ \vdots\\ y_n\\ \end{array} \right)\]

简记为 $X\boldsymbol{a}=\boldsymbol{y}$, 显然就有 $\boldsymbol{a}=X^{-1}\boldsymbol{y}$:

\[\left( \begin{array}{c} a_0\\ a_1\\ a_2\\ \vdots\\ a_n\\ \end{array} \right) =\left( \begin{matrix} 1& x_0& {x_0}^2& \cdots& {x_0}^n\\ 1& x_1& {x_1}^2& \cdots& {x_1}^n\\ 1& x_2& {x_2}^2& \cdots& {x_2}^n\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_n& {x_n}^2& \cdots& {x_n}^n\\ \end{matrix} \right) ^{-1}\left( \begin{array}{c} y_0\\ y_1\\ y_2\\ \vdots\\ y_n\\ \end{array} \right)\]

实际上这里的 $X$ 为范德蒙矩阵, 也就是说只要 $x_0,x_1,x_2,\ldots,x_n$ 互不相同, 那么 $X$ 就一定可逆.

2. 任务

我们的任务, 就是计算两个多项式的乘积: 对于

\[f\left( x \right) =a_0+a_1x+a_2x^2+\cdots +a_nx^n\]

和

\[g\left( x \right) =b_0+b_1x+b_2x^2+\cdots +b_mx^m\]

其乘积为:

\[h\left( x \right) =f\left( x \right) g\left( x \right) =c_0+c_1x+c_2x^2+\cdots +c_{n+m}x^{n+m}\]

我们希望能快速计算出系数 $c_0,c_1,…,c_{n+m}$.

3. 步骤

由于新系数有 $n+m+1$ 个, 所以我们将采样点也增加到 $n+m+1$ 个:

\[\left( \begin{matrix} 1& x_0& {x_0}^2& \cdots& {x_0}^n\\ 1& x_1& {x_1}^2& \cdots& {x_1}^n\\ 1& x_2& {x_2}^2& \cdots& {x_2}^n\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_{n+m}& x_{n+m}^2& \cdots& x_{n+m}^n\\ \end{matrix} \right) \left( \begin{array}{c} a_0\\ a_1\\ a_2\\ \vdots\\ a_n\\ \end{array} \right) =\left( \begin{array}{c} f_0\\ f_1\\ f_2\\ \vdots\\ f_{n+m}\\ \end{array} \right)\] \[\left( \begin{matrix} 1& x_0& {x_0}^2& \cdots& {x_0}^m\\ 1& x_1& {x_1}^2& \cdots& {x_1}^m\\ 1& x_2& {x_2}^2& \cdots& {x_2}^m\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_{n+m}& x_{n+m}^2& \cdots& x_{n+m}^m\\ \end{matrix} \right) \left( \begin{array}{c} b_0\\ b_1\\ b_2\\ \vdots\\ b_m\\ \end{array} \right) =\left( \begin{array}{c} g_0\\ g_1\\ g_2\\ \vdots\\ g_{n+m}\\ \end{array} \right)\] \[\left( \begin{matrix} 1& x_0& {x_0}^2& \cdots& {x_0}^{n+m}\\ 1& x_1& {x_1}^2& \cdots& {x_1}^{n+m}\\ 1& x_2& {x_2}^2& \cdots& {x_2}^{n+m}\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_{n+m}& x_{n+m}^{2}& \cdots& x_{n+m}^{n+m}\\ \end{matrix} \right) \left( \begin{array}{c} c_0\\ c_1\\ c_2\\ \vdots\\ c_{n+m}\\ \end{array} \right) =\left( \begin{array}{c} h_0\\ h_1\\ h_2\\ \vdots\\ h_{n+m}\\ \end{array} \right)\]

其中 $h(x_i)=f(x_i)g(x_i)$:

\[\left( \begin{array}{c} h_0\\ h_1\\ h_2\\ \vdots\\ h_{n+m}\\ \end{array} \right) =\left( \begin{array}{c} f_0\\ f_1\\ f_2\\ \vdots\\ f_{n+m}\\ \end{array} \right) \odot \left( \begin{array}{c} g_0\\ g_1\\ g_2\\ \vdots\\ g_{n+m}\\ \end{array} \right) =\left( \begin{array}{c} f_0g_0\\ f_1g_1\\ f_2g_2\\ \vdots\\ f_{n+m}g_{n+m}\\ \end{array} \right)\]

最终我们需要的系数向量为: $\boldsymbol{c}=X^{-1}\boldsymbol{h}$.

4. 快速傅里叶变换

但是如果用上面的方法, $X^{-1}$ 要算到猴年马月去了, 关键在于采样点的选取, 一方面需要使得 $X^{-1}$ 便于计算, 另一方面, 最好能让 $X\boldsymbol{a}$ 和 $X\boldsymbol{b}$ 的值能被快速计算出来, 有没有这样的好事呢？还真有。

我们先记 $n+m+1=N$, 那么采样点为 $x_i=\omega ^i$, 其中 $\omega$ 为单位1 的 $N$ 次根 $\mathrm{e}^{2\pi\mathrm{i}/{N}}$.

4.1. 关于 X

那么 $X$ 矩阵为

\[X=\left( \begin{matrix} 1& x_0& {x_0}^2& \cdots& {x_0}^{N-1}\\ 1& x_1& {x_1}^2& \cdots& {x_1}^{N-1}\\ 1& x_2& {x_2}^2& \cdots& {x_2}^{N-1}\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_{N-1}& x_{N-1}^{2}& \cdots& x_{N-1}^{N-1}\\ \end{matrix} \right) =\left( \begin{matrix} 1& 1& 1& \cdots& 1\\ 1& \omega& \omega ^2& \cdots& \omega ^{N-1}\\ 1& \omega ^2& \omega ^4& \cdots& \omega ^{2N-2}\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& \omega ^{N-1}& \omega ^{2N-2}& \cdots& \omega ^{\left( N-1 \right) ^2}\\ \end{matrix} \right)\]

其逆矩阵为

\[X^{-1}=\frac{1}{N}\left( \begin{matrix} 1& 1& 1& \cdots& 1\\ 1& \bar{\omega}& \bar{\omega}^2& \cdots& \bar{\omega}^{N-1}\\ 1& \bar{\omega}^2& \bar{\omega}^4& \cdots& \bar{\omega}^{2N-2}\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& \bar{\omega}^{N-1}& \bar{\omega}^{2N-2}& \cdots& \bar{\omega}^{\left( N-1 \right) ^2}\\ \end{matrix} \right)\]

这里的 $\bar{\omega}$ 为其共轭值。验证起来也很简单, 我们取 $X$ 的第 $i$ 行与 $X^{-1}$ 的第 $j$ 列相乘:

\[\begin{aligned} X\left[ i,: \right] X^{-1}\left[ :,j \right] &=\frac{1}{N}\left( \begin{matrix} 1& \omega ^i& \omega ^{2i}& \cdots& \omega ^{\left( N-1 \right) i}\\ \end{matrix} \right) \left( \begin{array}{c} 1\\ \bar{\omega}^j\\ \bar{\omega}^{2j}\\ \vdots\\ \bar{\omega}^{\left( N-1 \right) j}\\ \end{array} \right)\\ &=\dfrac{1+\omega ^{i-j}+\omega ^{2\left( i-j \right)}+\cdots +\omega ^{\left( N-1 \right) \left( i-j \right)}}{N}\\ &=\begin{cases} \dfrac{1-\omega ^{N(i-j)}}{1-\omega ^{i-j}}=0& i\ne j\\ 1& i=j\\ \end{cases}\\ \end{aligned}\]

4.2. 关于 Xa

上面的只是开胃菜, 我们如何快速计算 $X\boldsymbol{a}$ 呢？举个例子, 对于奇数阶多项式（N为偶数）:

\[f\left( x \right) =a_0+a_1x+a_2x^2+\cdots +a_{N-1}x^{N-1}\]

我们按奇偶进行分组:

\[\begin{aligned} f(x)&=\left( a_0+a_2x^2+\cdots +a_{N-2}x^{N-2} \right) +\left( a_1x+a_3x^3+\cdots +a_{N-1}x^{N-1} \right)\\ &=\left( a_0+a_2x^2+\cdots +a_{N-2}x^{N-2} \right) +x\left( a_1+a_3x^2+\cdots +a_{N-1}x^{N-2} \right)\\ \end{aligned}\]

然后构建新函数

\[G\left( x \right) =a_0+a_2x^1+\cdots +a_{N-2}x^{\frac{N}{2}-1}\] \[H\left( x \right) =a_1+a_3x^1+\cdots +a_{N-1}x^{\frac{N}{2}-1}\]

那么原来的 $f(x)$ 用新函数表示为

\[f(x)=G\left(x^{2}\right)+xH\left(x^{2}\right)\]

当我们带入 $x={\omega _N}^k\,\,( k<\frac{N}{2})$ 可得:

\[\begin{aligned} f({\omega_N}^k)&=G\left( {\omega _N}^{2k} \right) +{\omega_N}^kH\left( {\omega _N}^{2k} \right)\\ &=G\left( {\omega_{N/2}}^k \right) +{\omega _N}^kH\left( {\omega_{N/2}}^k \right)\\ \end{aligned}\]

而对于 $x={\omega _N}^k\,\,( k>\frac{N}{2})$ , 有

\[\begin{aligned} f({\omega_N}^k)&=G\left( {\omega _N}^{2k} \right) +{\omega_N}^kH\left( {\omega _N}^{2k} \right)\\ &=G\left( {\omega_N}^{2\left( k-N/2 \right)} \right) +{\omega _N}^kH\left( {\omega_N}^{2\left( k-N/2 \right)} \right)\\ &=G\left( {\omega _{N/2}}^{k-N/2} \right) -{\omega_N}^{k-N/2}H\left( {\omega _{N/2}}^{k-N/2} \right)\\ \end{aligned}\]

可以看到 $G\left( {\omega_{N/2}}^{k-N/2} \right)$ 和 ${\omega_N}^{k-N/2}H\left( {\omega_{N/2}}^{k-N/2} \right)$ 的值在上一步已经算过了, 而对于子问题 $G\left( {\omega_{N/2}}^{k-N/2} \right)$ 的计算, 我们可以继续使用上面的方法进行递归分治。

4.3. 代码实现

MATLAB 代码的实现如下：

function [y] = my_fft(a)
    n = size(a,2);
    if n == 1
        y = a;
        return
    end
    w =  complex(cos(2*pi/n),sin(2*pi/n));
    a_odd = a(1:2:end);
    a_even = a(2:2:end);
    g = my_fft(a_odd);
    h = my_fft(a_even);
    W = w.^(0:n/2-1);
    y = [g+W.*h, g-W.*h];
end

这里我们默认 $N$ 为 $2^k$，否则的话递归调用时会出现奇数长度的子问题，处理起来非常麻烦，最好的方式是在调用前补零，把数组长度扩充到2的幂次方。

4.4. 时间复杂度

对于一个长度为 $N$ 的数组，我们设 FFT 变换的计算量为 $T(N)$，那么有

\[\begin{aligned} T(N)&=N+2T(N/2)\\ &=N\log_2 N+ N T(1)\\ &= \mathcal{O}(N\log_2 N) \end{aligned}\]

4.5. 关于 FFT

为什么叫快速傅里叶变换呢？我们其实要计算的是 $\boldsymbol{a}$ 和 $\boldsymbol{b}$ 的卷积, 这里的左乘 $X$ 就相当于一个傅里叶变换, 这样的话, 我们就可以在傅里叶域进行点积, 然后再逆变换回来。

本来算卷积是 $\mathcal{O} ( N^2)$ 的复杂度, 但傅里叶变换的复杂度为 $\mathcal{O} \left( N\log N \right)$, 点积的复杂度为 $\mathcal{O} \left( N\right)$, 因此总复杂度为 $\mathcal{O} \left( N\log N \right)$, 节约了更多时间。

5. 流程总结

总结一下，如果我们快速计算两个多项式的乘积, 本质上是算其多项式系数的卷积 $\boldsymbol{a}\star \boldsymbol{b}$, 其流程如下:

  1. 先通过 FFT 将其转化到点值域: $\boldsymbol{f}=X\boldsymbol{a}$ 和 $\boldsymbol{g}=X\boldsymbol{b}$.
  2. 在点值域进行点积: $\boldsymbol{h}=\boldsymbol{f}\odot \boldsymbol{g}$.
  3. 再通过 FFT 的逆变换转化为系数表达形式: $\boldsymbol{c}=X^{-1}\boldsymbol{h}$.

6. 实验

我们可以用 MATLAB 中的 conv 函数验证我们的计算结果，或着也可以用内建的 fft 函数，但需要注意一点：MATLAB 内建 fft 函数中的 $w=\mathrm{e}^{-2\pi\mathrm{i}/{N}}$，因此算出来的结果会有一个共轭的区别，但对整体的结果没有影响。

另外需要注意的一点是：我们的 my_fft函数只能处理长度为2的幂次方的情况，所以要记得补零。

clc,clear
n = 7;
m = 5;
a = randi(10,1,n);
b = randi(10,1,m);

N = 2^ceil(log2(n+m));
a_ = [a, zeros(1,N-n)];
b_ = [b, zeros(1,N-m)];
A = my_fft(a_);
B = my_fft(b_);
C = A.*B;
c_ = real(my_ifft(C)/N);
c = c_(1:m+n-1)
c_conv = conv(a,b)

% A1 = fft(a_);
% B1 = fft(b_);
% C1 = A1.*B1;
% c_1 = real(ifft(C1))


function [y] = my_fft(a)
    n = size(a,2);
    if n == 1
        y = a;
        return
    end
    w = complex(cos(2*pi/n),sin(2*pi/n));
    a_odd = a(1:2:end);
    a_even = a(2:2:end);
    g = my_fft(a_odd);
    h = my_fft(a_even);
    W = w.^(0:n/2-1);
    y = [g+W.*h, g-W.*h];
end

function [y] = my_ifft(a)
    n = size(a,2);
    if n == 1
        y = a;
        return
    end
    w = complex(cos(2*pi/n),-sin(2*pi/n));
    a_odd = a(1:2:end);
    a_even = a(2:2:end);
    g = my_ifft(a_odd);
    h = my_ifft(a_even);
    W = w.^(0:n/2-1);
    y = [g+W.*h, g-W.*h];
end

7. 参考

  1. 快速傅里叶变换（FFT）超详解
  2. 快速傅里叶变换 - OI Wiki