有趣的算法03--多项式乘法

今天这篇文章主要讨论多项式乘法, 这大概是目前为止我写过的最为硬核的文章, 但是我会尽量用最为通俗的语言去解释, 相信各位有高中学历的同学们都能看懂.

多项式乘法在日常生活中的投资理财、市场预测、建筑施工等方面都有应用:

• 计算复利：在计算复利时，如果投资期限较长，涉及多次利率计算，就可以用多项式乘法的思路来理解。假设本金为 \(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\) 为各项系数。通过多项式乘法运算，可以根据当前的时间点预测未来一段时间的销售额，帮助商家合理安排库存、制定营销策略等。

同时, 两个整数的乘法也可以通过多项式乘法来求出. 由此可见, 多项式乘法在生活中是极其重要的.

问题描述

对于两个多项式:

希望求出多项式 \(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\) 上的值来唯一确定这个多项式的系数.

首先, 我先介绍一些后续需要的定义和定理:

1. 矩阵乘法: 若矩阵 \(A_{m \times n}\) 的第 \(i\) 行第 \(j\) 列元素记为 \(a_{ij}\) , 同时 \(B_{n \times s}[j, k] = b_{jk}\) , 则

\[(A \cdot B)[i, j] = \sum_{k = 1}^n a_{ik}b_{kj}\]

1. 矩阵求逆: 一个矩阵的逆可以理解为一个矩阵的”倒数”: 一个矩阵(数)乘以它的逆(倒数)等于单位矩阵( \(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)

首先, 让我们回忆一些后续需要的高中知识:

1. 复平面: 复平面中每个点都对应一个复数, 其中 \((a, b)\) 对应 \(a + bi\). 值得注意的是, 我们也可以使用极坐标表示一个点: \((\rho, \theta)\) 表示复数 \(\rho e^{\theta i}\).
2. 复数乘法: 极坐标中 \((\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\) 有以下性质:

1. \[\omega^{i} = \omega^{i + n}\]
2. \(\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\) 处的取值:

因此, 欲求 \(P(\omega^{0}), P(\omega^{1}), \dots, P(\omega^{n-1})\) , 只需求

和

而这是 \(2\) 个 \(\frac{n}{2}\) 阶的多项式在 \(\frac{n}{2}\) 个点上的值( \(\omega^{i} = \omega^{n+i}\) ), 因此可以得出时间复杂度 \(T(n)\) 的递推公式

好耶! 我们的多项式乘法终于变成 \(O(n\log n)\) 的算法啦!