一元多项式

多项式的形式很简单，但其功能十分强大，在数学上和工程上有着极其重要的地位。数学上任何闭区间上的连续函数，都可以用多项式逼近到任意精度。在 CAD、动画制作等工程领域中，都是在用各种多项式来描述复杂的曲线曲面。在控制系统中，人们也习惯使用多项式来描述系统的传递函数，分析系统的稳定性。

一个一元 $n$ 次多项式 $P(x)$，可以写成如下的形式:

$$ \begin{equation}\label{f1} P(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0 \end{equation} $$

其中，$a_i, i \in \{0,1,\cdots,n\}$ 是多项式的系数，特别的 $a_n \neq 0$。$n$ 就是多项式 $P(x)$ 的次数，$x$ 是其自变量。因为只有一个自变量 $x$，所以这里强调是一元的。多项式也可以是多元的，即有多个自变量，以后我们可能会有专题来讨论它。这里我们只关注实系数的一元多项式。

1. 求值

我们先来研究一下如何高效的计算多项式的函数值。如果我们完全按照式($\ref{f1}$)的形式逐项计算 $a_i x^i$ 再求和，那么我们需要计算 $(n + 1)^2/2$ 次乘法和 $n$ 次加法。即使复用 $x^i$ 的中间结果，也需要 $2n-1$ 次乘法和 $n$ 次加法。实际上，如果我们对式($\ref{f1}$)做如下的变形，就可以只用 $n$ 次乘法和 $n$ 次加法完成计算。人们常把这种计算方法称为嵌套乘法或Horner 方法。

$$ \begin{equation}\label{f2} P(x) = \left(\cdots\left(a_n x + a_{n-1}\right)x + \cdots + a_1\right)x + a_0 \end{equation} $$

式($\ref{f2}$)相比于式($\ref{f1}$)而言，在数值计算上也更稳定。对于式($\ref{f1}$)，如果多项式的次数很高而且自变量 $x$ 远大于 1，那么我们在计算 $x^n$ 的时候就很有可能溢出。如果 $x$ 很小十分接近 0，那么 $x^n$ 就很可能因为浮点数的舍入误差变得毫无意义。而式$\ref{f2}$的计算过程就要安全很多，每次乘法计算之后都会附带一次加法运算，是有机会调整中间数值的，可以有效避免过大或过小的乘法运算。相关实现参考Polynomial::Evaluate。

2. 求根

求一个多项式的根，即方程 $a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0 = 0$，在工程上是一种基本需求。比如通过控制系统传递函数的零极点来分析系统的稳定性，其中零极点就是传递函数分子和分母上多项式方程的根。我们在第一部分中介绍的迭代方法一次只能收敛到一个根上。而一般的一个 $n$ 次多项式可以有 $n$ 个根(可能有重根或复数根)。人们研究出了一种称为收缩(deflation)的约化技术，可以解决这个问题，利用逐一求出所有的根。

对于多项式 $P(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0$，我们令:

$$ \begin{equation}\label{f3} \begin{cases} b_n = a_n & \\ b_k = a_k + b_{k+1}x_0, & k = n-1, n-2, \cdots, 1, 0 \end{cases} \end{equation} $$

有多项式 $Q(x) = b_nx^{n-1} + b_{n-1}x^{n-2} + \cdots + b_2 x + b_1$，使得:

$$ \begin{equation}\label{f4} P(x) = (x - x_0)Q(x) + b_0 \end{equation} $$

上述定理很容易证明，如果我们将式($\ref{f4}$)展开有:

$$ \begin{aligned} (x - x_0)Q(x) + b_0 & = (x - x_0)(b_nx^{n-1} + b_{n-1}x^{n-2} + \cdots + b_2 x + b_1) + b_0 \\ & = (b_n x^n + b_{n-1}x^{n-1} + \cdots + b_2x^2 + b_1x) \\ & \quad -(b_n x_0 x^{n-1} + \cdots + b_2 x_0 x + b_1 x_0) + b_0 \\ & = b_n x^n + (b_{n-1} - b_nx_0)x^{n-1} + \cdots + (b_0 - b_1x_0) \end{aligned} $$

我们将式($\ref{f3}$)稍作变形，有 $b_n = a_n$, $b_k - b_{k+1}x_0 = a_k$，代入上式有$(x - x_0)Q(x) + b_0 = P(x)$。□

上述定理有两个重要推论 $P(x_0) = b_0$ 和 $P'(x_0) = Q(x_0)$。其中 $P(x_0) = b_0$ 很显然只需将 $x_0$ 代入式($\ref{f4}$) 就能得到。根据求导法则有:

$$ P'(x) = Q(x) + (x - x_0)Q'(x) $$

将 $x_0$ 代入上式，就有 $P'(x_0) = Q(x_0)$。根据这两个推论和式($\ref{f3}$)，我们可以通过嵌套乘法同时求出 $P(x_0), P'(x_0)$，相关实现参考Polynomial::Horner。然后套用牛顿法，如果收敛的话，似乎就能找到一个根。如果 $x_0$ 就是 $P(x) = 0$ 的一个根，那么有 $b_0 = 0$，剩下的根好像都在方程 $Q(x) = 0$ 中。对 $Q(x)$ 再套用收缩约化，如此重复下去，我们可能就能找到多项式所有的根。

看起来我们找到了一种求解多项式所有根的算法。很不幸这个“算法”有很多缺陷。其一，我们知道在一些情况下多项式是存在复数根的，该“算法”的整个迭代过程都是在实数空间下的运算，遇到复数根就挂了。其二，该“该算法”在数值上会很不稳定，因为迭代求根是存在数值误差的，所以 $Q(x)$ 的各项系数并不是那么精确，会带来累积误差，导致最后求出的若干个根误差很大甚至是错误的结果。

我们会在下一节中介绍可以处理复数的 Müller 方法以及一些其它更稳定的求根方法。下面我们来研究一个特别的多项式。

3. 一元二次多项式

在中学阶段我们就已经背过，一个一元二次多项式方程 $ax^2 + bx + c = 0, a \neq 0$ 的求根公式:

$$ \begin{equation}\label{f5} x_1 = \frac{-b + \sqrt{b^2 - 4ac}}{2a} \qquad x_2 = \frac{-b - \sqrt{b^2 - 4ac}}{2a} \end{equation} $$

现在我们来看一个具体的方程 $(x - 10^4)(x - 10^{-4}) = 0$，显然它有两个根 $x_1 = 10000, x_2 = 0.0001$。展开有 $x^2 - 10000.0001 x + 1 = 0$。按照如下左侧的 python3 脚本计算一遍，会发现 x1 正好是我们要的结果，x2 稍微有一些误差。

        import math
        a = 1
        b = -10000.0001
        c = 1

        x1 = (-b + math.sqrt(b*b - 4*a*c)) / (2 * a)
        x2 = (-b - math.sqrt(b*b - 4*a*c)) / (2 * a)
        print(f"x1 = {x1}, x2 = {x2}")
        ----------------------------------------------
        x1 = 10000.0, x2 = 0.00010000000020227162

        abs_e1 = abs(x1 - 10000)
        rel_e1 = abs_e1 / 10000
        print(f"x1 绝对误差: {abs_e1}, 相对误差: {rel_e1}")
        
        abs_e2 = abs(x2 - 0.0001)
        rel_e2 = abs_e2 / 0.0001 
        print(f"x2 绝对误差: {abs_e2}, 相对误差: {rel_e2}")
        ----------------------------------------------
        x1 绝对误差: 0.0, 相对误差: 0.0
        x2 绝对误差: 2.0227161688212564e-13, 相对误差: 2.0227161688212564e-09

在这个例子中，我们故意构造了两个绝对值比较悬殊的根。如此 $b^2$ 就远大于 $4ac$ 这导致 $\sqrt{b^2 - 4ac} \approx |b|$，所以在计算 $x_2$ 的时候，分子上是两个绝对值接近的数相减，这就要面对浮点数舍入误差。而 $x_1$ 则是两个大数相加，舍入误差的影响就相对要小一些。 python 中浮点数都是双精度的，也就是 C/C++ 中的 double 类型，x2 的相对误差就已经达到了 $10-9$ 的量级。如果是 float 类型的，这点误差就可能是个错误了。

对于 $x_2$ 我还是可以通过分子有理化的方式，得到一个数值上更精确的解:

$$ x_2 = \frac{-b - \sqrt{b^2 - 4ac}}{2a} \cdot \frac{-b + \sqrt{b^2 - 4ac}}{-b + \sqrt{b^2 - 4ac}} = \frac{b^2 - (b^2 - 4ac)}{2a(-b + \sqrt{b^2 - 4ac})} = \frac{-2c}{b - \sqrt{b^2 - 4ac}} $$

如此，我们将在分子上两个相近数的减法，变成了分母上两个大数的加法。这样舍入误差的影响就小了。对称的，我们也可以对 $x_1$ 做分子有理化，得到如下求根公式。

$$ \begin{equation}\label{f6} x = \frac{-2c}{b \pm \sqrt{b^2 - 4ac}} \end{equation} $$

结果上也很对称，$x_1$ 原本是分子上的大数加法，变成了分母上的减法，可以预见按照上式计算 $x_1$ 就要面对舍入误差了。正如下面的代码片段和日志所示。所以实际应用求根公式的时候，需要根据 $b$ 的符号在式($\ref{f5}$)和式($\ref{f6}$)中二选一。

        import math
        a = 1
        b = -10000.0001
        c = 1

        x1 = (-2 * c) / (b + math.sqrt(b*b - 4*a*c))
        x2 = (-2 * c) / (b - math.sqrt(b*b - 4*a*c))
        print(f"x1 = {x1}, x2 = {x2}")
        ----------------------------------------------
        x1 = 9999.999979772838, x2 = 0.0001

        abs_e1 = abs(x1 - 10000)
        rel_e1 = abs_e1 / 10000
        print(f"x1 绝对误差: {abs_e1}, 相对误差: {rel_e1}")
        
        abs_e2 = abs(x2 - 0.0001)
        rel_e2 = abs_e2 / 0.0001 
        print(f"x2 绝对误差: {abs_e2}, 相对误差: {rel_e2}")
        ----------------------------------------------
        x1 绝对误差: 2.0227162167429924e-05, 相对误差: 2.0227162167429924e-09
        x2 绝对误差: 0.0, 相对误差: 0.0

我们在这里讨论一元二次方程的主要目的是，强调不同的计算姿势下，舍入误差的影响是不一样的。在设计或者选型算法的时候，不能只看数学形式上是否正确，还要仔细考虑一下大数和小数的问题。因为舍入误差导致的问题往往很难定位。

4. 完

本节我们介绍了用于高效计算多项式函数值的嵌套乘法。然后讨论了利用收缩约化技术求解多项式方程所有根的方法，真正可以应用到工程中的方法将放到下一节介绍。最后，我们讨论了一元二次多项式方程的求根公式，用一个具体的例子分析了它的舍入误差。

本节中我更想强调的是，嵌套乘法关于计算溢出的考虑，以及最后一元二次方程求根公式中关于舍入误差的考虑。差之毫厘谬以千里，在数值计算中我们需要时刻对这两类问题保持警惕。