今天这篇文章主要讨论多项式乘法, 这大概是目前为止我写过的最为硬核的文章, 但是我会尽量用最为通俗的语言去解释, 相信各位有高中学历的同学们都能看懂.
多项式乘法在日常生活中的投资理财、市场预测、建筑施工等方面都有应用:
- 计算复利：在计算复利时，如果投资期限较长，涉及多次利率计算，就可以用多项式乘法的思路来理解。假设本金为 $P$ ，年利率为 $r$ ，投资 $n$ 年，每年的本息和计算就类似于多项式展开。第一年本息和为 $P(1 + r)$ ，第二年为 $P(1 + r)(1 + r)=P(1 + r)^2$ ，以此类推，第 $n$ 年为 $P(1 + r)^n$ 。这里 $(1 + r)^n$ 展开后就是一个关于 $r$ 的 $n$ 次多项式，通过多项式乘法的原理可以清晰地看到每年利息对最终收益的影响，帮助投资者准确计算收益。
- 销售趋势分析：商家在分析某种商品的销售趋势时，可能会根据过去几个月或几年的销售数据来建立一个多项式模型。例如，以时间 $t$ 为变量，销售额 $S$ 可能可以表示为 $S = a_0 + a_1t + a_2t^2+\cdots+a_nt^n$ 的形式，其中 $a_i$ 为各项系数。通过多项式乘法运算，可以根据当前的时间点预测未来一段时间的销售额，帮助商家合理安排库存、制定营销策略等。
同时, 两个整数的乘法也可以通过多项式乘法来求出. 由此可见, 多项式乘法在生活中是极其重要的.
问题描述
对于两个多项式:
\[\begin{align}& A(x) = a_0 + a_1x^1 + \dots + a_nx^n \\ & B(x) = b_0 + b_1x^1 + \dots + b_mx^m\end{align}\]希望求出多项式 $C(x) = A(x) \times B(x)$ 的各项系数.
我们这里不妨假设两个输入多项式最多都是 $n$ 阶的, 而输出的多项式最多是 $2\cdot n$ 阶的. 由于不同多项式之间的区别主要体现在系数上, 因此我们可以使用一个 $n$ 阶的系数向量来表示一个多项式.
A[i] = a_i => A(x) = A[0] + A[1] * x ** 1 + ...
注意: 这里我们认为乘法的时间复杂度是 $O(1)$ .
算法I: 循环
对于这样一个问题, 最朴素的方式就是通过两个循环分别计算 $C(x)$ 中各项系数. 伪代码如下:
def poly_mul(A, B):
C[2n] = {0}
for i in range(n):
for j in range(n):
C[i + j] += A[i] * B[j]
return C
显而易见的是, 这个算法的时间复杂度是 $O(n^2)$ .
算法II: 多项式插值法
下面, 我们先介绍一个需要使用的定理:
定理: 对于任何一个 $n$ 阶多项式, 我们可以通过确定这个多项式在 $n + 1$ 以上不同的 $x$ 上的值来唯一确定这个多项式的系数.
首先, 我先介绍一些后续需要的定义和定理:
\[(A \cdot B)[i, j] = \sum_{k = 1}^n a_{ik}b_{kj}\]
- 矩阵乘法: 若矩阵 $A_{m \times n}$ 的第 $i$ 行第 $j$ 列元素记为 $a_{ij}$ , 同时 $B_{n \times s}[j, k] = b_{jk}$ , 则
- 矩阵求逆: 一个矩阵的逆可以理解为一个矩阵的”倒数”: 一个矩阵(数)乘以它的逆(倒数)等于单位矩阵( $1$ ). 值得注意的是, 只有行列式不为 $0$ 的方阵才可以求逆.
证明如下:
\[\begin{align*}&\begin{cases}a_0 + a_1x_1 + \dots + a_nx_1^n = y_1\\a_0 + a_1x_2 + \dots + a_nx_2^n = y_2\\ \dots\\ a_0 + a_1x_{n+1} + \dots + a_nx_{n+1}^n = y_{n+1}\end{cases} \implies \begin{pmatrix}1 & x_1 & \dots & x_1^n \\ 1 & x_2 & \dots & x_2^n \\ &\dots \\ 1 & x_{n+1} & \dots &x_{n+1}^n\end{pmatrix}\cdot\begin{pmatrix}a_0\\a_1\\\dots\\a_n\end{pmatrix} = \begin{pmatrix}y_1\\y_2\\\dots\\y_{n+1}\end{pmatrix} \\ &A = \begin{pmatrix}1 & x_1 & \dots & x_1^n \\ 1 & x_2 & \dots & x_2^n \\ &\dots \\ 1 & x_{n+1} & \dots &x_{n+1}^n\end{pmatrix}\text{的行列式不为 }0\text{ , 因此可以求 }A^{-1}\text{ . 我们称形如 }A\text{ 的矩阵为范德蒙德矩阵} \\ &\therefore \begin{pmatrix}a_0\\a_1\\\dots\\a_n\end{pmatrix} = A^{-1}\cdot \begin{pmatrix}y_1\\y_2\\\dots\\y_{n+1}\end{pmatrix}\end{align*}\]
这个定理是比较直观的, 实质上, 我们所说的”两点确定一条直线”, “三点确定一个平面“都是这个定理在低维上的体现.
那么, 对于多项式 $C = A \times B$ , 我们就可以分别确定 $A$ 和 $B$ 在 $2n + 1$ 个不同 $x$ 上的取值 $A(x_1), A(x_2), \dots, A(x_{2n+1})$ 和 $B(x_1), B(x_2), \dots, B(x_{2n+1})$ , 并将它们对应相乘, 就可以得出 $C(x_1), C(x_2), \dots, C(x_{2n+1})$, 因此我们就可以确定 $C$ 的系数啦!
求 $C$ 的系数向量也很简单, 只需要按照证明, 写出任意 $2n + 1$ 个 $x$ , 并求出它们在 $C$ 上的值向量 $Y$. 接着, 我们求出这些 $x$ 对应的那个 $2n+1$ 阶的范德蒙德矩阵 $A$ , 并求它的逆 $A^{-1}$ . 最后, 求出 $C$ 的系数向量 $C = A^{-1}\cdot Y$. 伪代码如下:
def poly_mul(A, B):
X = [i for i in range(2 * n + 1)]
Y = [0 for i in range(2 * n + 1)]
for i in range(2 * n + 1): # O(n ^ 2)
a, b = 0, 0
for j in range(n):
a += A[j] * X[i] ** j
b += B[j] * X[i] ** j
Y[i] = a * b
V = [[0 for _ in range(2 * n + 1)] for _ in range(2 * n + 1)]
for i in range(2 * n + 1): # O(n ^ 2)
for j in range(2 * n + 1):
V[i][j] = X[i] ** j
C = get_inverse(V) * Y
return C
由伪代码可以看出, 改算法至少是 $O(n^2)$ 的时间复杂度.
实际上,
get_inverse(V)的时间复杂度是 $O(n^3)$ . 因此, 这个算法的时间复杂度甚至不如算法I.
真是太糟糕了! 难道这么巧妙的方法就没有用武之地嘛! 要是我们能够快速找到一个范德蒙德矩阵和它的逆就好了.
算法III: 快速傅里叶变换(FFT)
首先, 让我们回忆一些后续需要的高中知识:
- 复平面: 复平面中每个点都对应一个复数, 其中 $(a, b)$ 对应 $a + bi$. 值得注意的是, 我们也可以使用极坐标表示一个点: $(\rho, \theta)$ 表示复数 $\rho e^{\theta i}$.
- 复数乘法: 极坐标中 $(\rho_1, \theta_1)$ 和 $(\rho_2, \theta_2)$ 表示的复数乘积为 $(\rho_1\cdot\rho_2, \theta_1+\theta_2)$. 这个定理很好证明, 可以使用三角函数和角公式推出, 这里就先不赘述了.
让我们考虑一个复数 $\omega = e^{2\pi i / n}$. 这个复数在极坐标中对应的点是 $(1, 2\pi /n)$ . 同时, $\omega^2$ 对应的点是 $(1, 4\pi /n)$, $\omega ^ k$ 对应的点是 $(1, 2k\pi / n)$ . 可以发现, $1, \omega, \omega^2, \dots, \omega^{n-1}$ 对应的是复平面中单位圆上从x轴开始不断逆时针旋转 $2\pi /n$的一系列点.
显然, $\omega$ 有以下性质:
- $\omega^{i} = \omega^{i + n}$
- $\bar{\omega^i} = \omega^{n-i}$, 其中 $\bar{a}$ 是 $a$ 的共轭复数.
同时, 考虑范德蒙德矩阵 $V_{ij} = (\omega^i)^j = \omega^{ij}$, 经过一系列精彩的特征值分解可以得出: $V^{-1} = \frac{1}{n}\bar{V}$ .
哇! 这样我们就能快速得到一个范德蒙德行列式和它的逆啦! 考虑算法II, 还有一个问题: 算出 $C(\omega), C(\omega^2), \dots$ 仍然是 $O(n^2)$ 的算法.
怎么办! 难道我们的努力又又付之东流了吗!!!
当然不是啦! 让我们考虑一个多项式 $P$ 在 $x$ 处的取值:
\[\begin{align}P(x) = & a_{0} + a_{1}x + a_{2}x^{2} + \dots \\= & (a_{0} + a_{2}x^{2} + a_{4}(x^{2})^{2} + \dots) + (a_{1} + a_{3}x^{2} + a_{5}(x^{2})^{2} + \dots)x\\= & P_{even}(x^{2}) + P_{odd}(x^{2})x\end{align}\]因此, 欲求 $P(\omega^{0}), P(\omega^{1}), \dots, P(\omega^{n-1})$ , 只需求
\[P_{even}(\omega^{0}), P_{even}(\omega^{2}), \dots, P_{even}(\omega^{2n-2})\]和
\[P_{odd}(\omega^{0}), P_{odd}(\omega^{2}), \dots, P_{odd}(\omega^{2n-2})\]而这是 $2$ 个 $\frac{n}{2}$ 阶的多项式在 $\frac{n}{2}$ 个点上的值( $\omega^{i} = \omega^{n+i}$ ), 因此可以得出时间复杂度 $T(n)$ 的递推公式
\[T(n) = 2T\left( \frac{n}{2} \right) + \Theta(n)\implies T(n) = O(n\log n)\]好耶! 我们的多项式乘法终于变成 $O(n\log n)$ 的算法啦!