多项式算法 C++

快速稳定的多项式求根器 - 第一部分

Henrik Vestermark

5.00/5 (10投票s)

2023年10月10日

CPOL

17分钟阅读

14339

使用牛顿方法实现快速、鲁棒且可靠的多项式求根器的实践

引言

牛顿法求多项式根是最流行和简单的方法之一。牛顿方法使用以下算法逐步找到越来越接近根的值。

一次找到一个根。存在其他方法可以同时逐步逼近所有根。这些方法如 Aberth-Ehrlich 或 Durand-Kerner。然而，它们存在其他问题，使得实现起来不太理想。当然，还有许多其他方法可以考虑。其中有 Halley、Householder 3^rd order、Ostrowski、Laguerre、Graeffe、Jenkins-Traub（很可能最著名）、特征值方法等等。所有这些方法都有快速稳定的版本。读者可以参考 [1]，其中介绍了 20 多种不同的方法及其实现。现在，我将重点介绍使用牛顿方法实现鲁棒稳定求根器的实践。我们还将要求多项式具有复数系数。算法与多项式具有实数（第二部分）还是复数系数（本文）无关。

当前任务

使用牛顿法寻找多项式根通常易于实现

通常，您会经历以下步骤。

消除根为零的简单根。
设置牛顿迭代
1. 找到一个合适的初始猜测值
2. 使用 Horner 方法评估 P(x_n) 和 P’(x_n)
3. 找到步长 P(x_n)/P’(x_n)
4. 计算下一个 x_n+1
5. 重复 b-d 直到满足停止条件
6. 将新找到的根在 P(x) 中进行除法，以计算新的归约多项式 P_new(x)
7. 重复 a-f 直到剩下一次或二次多项式
直接求解一次或二次多项式

对于许多其他求根方法，您会经历或多或少相同的步骤。

牛顿方法的缺陷

牛顿方法本身并不一定稳定，但需要额外的代码来处理在实现鲁棒、快速、稳定解决方案时遇到的经典陷阱。仅看上面的公式（2），很明显当 P’(x) 为零或接近零时我们会遇到问题。但这并不是你会遇到的唯一问题。有时，如果我们遇到局部最小值，步长 P(x)/P’(x) 可能会将搜索引向远离任何根的位置。从哪里开始搜索根也很重要，因为上述算法只有在您离根足够近时才会收敛到根。牛顿方法的收敛阶数为 2，这意味着每次迭代会使正确的位数翻倍。但当尝试寻找包含多个根的解时，例如 (x-2)²(x-3)=0，收敛速率会降至线性速率，需要更多迭代才能找到根。

我们需要解决多重根问题，并确保在这些情况下也能保持 2 阶收敛。在某个时候，我们需要确定何时停止搜索并满意于精度。我们的目标是将 IEEE754 浮点标准的能力发挥到极致，如同 C++ double 类型所实现的。如果我们放宽停止条件，不准确之处将传播到其他根，使它们离真实根越来越远。最后，当找到一个根后，我们需要将该根从多项式中除掉，并对新的归约多项式重复搜索。进行综合除法时，您可以在前向除法、后向除法或复合除法之间进行选择。选择可能会影响根的精度。

多重根问题

考虑多项式

P(x)=(x-1)(x-2)(x-3)(x-4)=x⁴-10x³+35x²-50x+24

如下所示，根是分开的。

图 1. 分开的根

使用起始点 0.5，牛顿迭代如下逼近第一个根

	x	P(x)
初始猜测	0.5
1	0.798295454545455	6.6E+00
2	0.950817599863883	1.7E+00
3	0.996063283034122	3.2E-01
4	0.999971872651986	2.4E-02
5	0.999999998549667	1.7E-04
6	0.999999999999999	8.7E-09
7	1.000000000000000	7.1E-15

正如我们所见，经过 7 次迭代就得到了第一个根 x=1。我们也注意到，在第二次迭代 x₂=0.95 之后，我们每次迭代都大致使第一个根的正确位数翻倍。每迭代使正确位数翻倍的迭代方法被称为具有 2 阶收敛。

现在我们更改多项式，并在 x=1 处引入一个重根

P(x)=(x-1)² (x-3)(x-4)=x⁴-9x³+27x²-31x+12

使用相同的起始点 x=0.5，收敛速度慢得多，经过 27 次迭代后，对于第一个根 x=1 没有任何改进，结果仅对前 8 位数字近似准确。

	x	P(x)
初始猜测	0.5
1	0.713414634146341	2.2E+00
2	0.842942878437970	6.2E-01
3	0.916937117337937	1.7E-01
…
10	0.999306565270595	1.2E-05
…
20	0.999999322514237	1.1E-11
…	0.999999661405383	2.8E-12
27	0.999999996306426	0.0E+00

这里到底发生了什么？

如果 P(x)=(x-1)² (x-3)(x-4)，那么 P'(x)=(x-1)(4x²-23x+31)

根 x=1 既是原始多项式 P(x) 的根，也是 P’(x) 的根。请看下图。

图 2. x=1 处的重根

在牛顿迭代中，当 P(x) 和 P’(x) 都趋向于 0 时，会在计算牛顿迭代的下一个 x_n+1 的精度中引入舍入误差。为说明起见，我们重复迭代步骤，但包括 P’(x)。此外，我们还引入了收敛速率 q。

	x	P(x)	P'(x)	q
初始猜测	0.5
1	0.713414634146341	2.2E+00	-1.0E+01
2	0.842942878437970	6.2E-01	-4.8E+00	1.3
3	0.916937117337937	1.7E-01	-2.3E+00	1.2
…
10	0.999306565270595	1.2E-05	-1.7E-02	1.1
…
20	0.999999322514237	1.1E-11	-1.6E-05	1.0
…
25	0.999999976999021	1.1E-14	-5.2E-07	1.0
26	0.999999996306426	5.3E-15	-2.8E-07	-
27	0.999999996306426	0.0E+00	-4.4E-08	-

我们注意到几点；收敛速率 q 比第一个例子慢得多；约 2 与约 1.1。此外，我们可以看到每次迭代，根的收敛都以线性因子 2 进行，而不是我们从第一个例子中的二次因子 2 所预期的。

对于更高阶的重根，情况会更糟。例如

如果 P(x)=(x-1)³ (x-4)，那么 P'(x)=(x-1)² (4x-13)

经过 31 次迭代，我们得到 x=0.999998662746209，仅对前 5 位数字近似准确。

	x	P(x)	P'(x)	q
初始猜测	0.5
1	0.659090909090909	4.4E-01	-2.8E+00
2	0.768989234449761	1.3E-01	-1.2E+00	1.2
…
29	0.999995827925540	4.4E-16	-2.9E-10	1.0
30	0.999998662746209	4.4E-16	-1.6E-10	-
31	0.999998662746209	0.0E+00	-1.6E-11	-

在三重重根 x=1 附近，P'(x) 在多重根周围的相当宽的半径内非常平坦，接近 0。

图 3. x=1 处的三重根

如何处理牛顿迭代中的多重根问题？

为了克服牛顿步长减小的这个问题，我们需要将其乘以因子 m，因此我们使用了修改后的牛顿迭代。

其中 m 是根的重数。这些都是众所周知的内容。挑战在于如何在现实生活中找到 m。

求根器的合适起始点

为了使迭代方法更快地收敛到多项式根，重要的是我们以某种方式从一个位于根附近的合适点开始。幸运的是，许多人已经研究过这个领域，并且 J.McNamee 在 [8]《Numerical Methods for Roots of Polynomials》中概述了 45 种以上创建根先验界限的方法。大多数先验界限是用来找到所有根所在的圆的半径。少数也涉及最小模根所在的圆的半径。这对于一次找到一个根的方法非常有用，其策略是按模大小递增的顺序找到根。

最小模根的先验界限

我所见的大多数求根器实现都不过多关注起始点。例如 [6] Grant-Hitchins 使用固定的起始点 (0.001+i0.1)。而不是固定的起始点，我更推荐 Madsen [2] 中实现的起始点。在该方法中，我们找到起始点 z₀，使得最小模根位于该圆的外部

最小根位于复平面上半径为 |z₀| 的圆的外部。

考虑多项式

P(x)=(x-1)(x+2)(x-3)(x-4)= x⁴+2x³-13x²-14x+24

上述公式给出了一个起始点 z₀=0.68，接近最近的根 x=1。

现在考虑多项式

P(x)= (x-1)(x+1000)(x-2000)(x+3000)= x⁴+1999x³-5002E3x²-5995e6x+6E9

上述公式给出了 z₀= 0.5（最近的根 x=1）。

找到第一个根 x=1 后，归约后的多项式为 P(x)= (x+1000)(x-2000)(x+3000)=x³+2E3x²-5E6x-6E9，上述公式给出了一个新起始点，用于搜索归约后的多项式，为 z₀=600（最近的根 x=1,000）。

该算法计算了合理且合适的起始点用于我们的根搜索。

    // Compute the next starting point based on the polynomial coefficients
    // A root will always be outside the circle from the origin and radius min
    auto startpoint = [&](const vector<complex<double>>& a)
    {
    const size_t n = a.size() - 1;
    double a0 = log(abs(a.back()));
    double min = exp((a0 - log(abs(a.front()))) / static_cast<double>(n));

    for (size_t i = 1; i < n; i++)
        if (a[i] != complexzero)
        {
            double tmp = exp((a0 - log(abs(a[i]))) / static_cast<double>(n - i));
            if (tmp < min)
                min = tmp;
        }

    return min*0.5;
    };

复数点上的多项式求值

大多数求根方法都需要我们在某一点上求多项式的值。

求多项式 P(z) 的值，其中

P(z)=a_n zⁿ+a_n-1z^n-1,…,a₁ z+a₀

我们使用 Horner [4] 方法，其递推关系为

b_n=a_n
b_k=b_k-1z+a_k, k=n-1,…,0
P(z)=b₀

此递推关系的最后一项 b₀ 即为 P(z) 的值。Horner 方法一直被认为是求多项式在给定点的值的最有效方法。该算法适用于系数为实数或复数。

    // Evaluate a polynomial with complex coefficients a[] at a complex point z and
    // return the result
    // This is the Horner's methods
    auto horner = [](const vector<complex<double>>& a, const complex<double> z)
    {
        const size_t n = a.size() - 1;
        complex<double> fval=a.front();
        eval e;

        for (size_t i = 1; i <= n; i++)
            fval = fval * z + a[i];

        e = { z, fval,abs(fval) };
        return e;
    };

求根的合适停止准则

在 [8] 中，他们讨论了许多不同的计算合适的停止准则的技术。另请参阅 [1]。许多求根器可以使用 Adams [5] 或 Hitching [6] 的方法来为具有实系数或复系数的多项式找到合适的停止准则。

在进行迭代方法时，您最终需要考虑要为求根器应用何种停止准则。由于大多数迭代求根器都使用多项式求值来进展，因此只有当 P(z) 的求值足够接近 0 时，我们才接受该根并停止搜索，这是很自然的。

算术运算中的误差

J.H.Wilkinson 在 [7] 中表明，执行代数运算中的误差界限为

ε<½β^1-t, β 是基数，t 是精度（假设四舍五入到最近）

注意 ½β^1-t= β^-t

对于 Intel 微处理器系列和 IEE754 标准的浮点运算，β= 2，对于 64 位浮点算术或 2^-53，t=53。

一个简单的上限

使用上述信息，可以为 n 次多项式找到一个简单的上限。

	多项式
*操作次数*	实系数	复系数
实数点	\|a_o\|·2n·2^-53	\|a_o\|·4n·2^-53
复数点	\|a_o\|·4n·2^-53	\|a_o\|·6n·2^-53

一个更好的上限

在此类别中，我们有 Adams [5] 和 Grant & Hitchins [6] 等人的多项式停止准则。

多项式求根器通常可以处理具有实数或复数系数的多项式，这些多项式在实数或复数上求值。由于 Adams 的停止准则用于实系数多项式，我们将使用 Grant & Hitchins 的界限，它相似但适用于复系数多项式。

Grant & Hitchins 复系数多项式的停止准则

复系数多项式 z_n 在复数点 z 上使用 Horner 方法求值。Grant 和 Hitchins [6] 使用递推关系推导出多项式求值的误差上限，其中 z=x+iy，复系数表示为 a_n+ib_n。

c_n=a_n, d_n=b_n, g_n=1, h_n=1;
c_k=xc_k+1-yd_k+1+a_k, k=n-1,...,0
d_k=yc_k+1+xd_k+1+b_k
g_k=|x|(g_k+1+|c_k+1|)+|y|(h_k+1+|d_k+1|)+|a_k|+2|c_k|
h_k=|y|(g_k+1+|c_k+1|)+|x|(h_k+1+|d_k+1|)+|b_k|+2|d_k|

现在误差是 (g₀+ih₀)ɛ，其中 ɛ= ½ β^1-t。现在由于递推关系本身引入了误差 [6]，我们通过加上递推关系中舍入误差的上限来保护计算，因此我们得到复数点上复数多项式求值的界限

e=(g₀+ih₀ )ε(1+ε)⁵ⁿ，其中 ε=½β^1-t

还有其他有用的方法可以考虑，请参阅 [1]。

    // Calculate an upper bound for the rounding errors performed in a
    // polynomial with complex coefficient a[] at a complex point z.
    // (Grant & Hitchins test)
    auto upperbound = [](const vector<complex<double>>& a, complex<double> z)
    {
        const size_t n = a.size() - 1;
        double nc, oc, nd, od, ng, og, nh, oh, t, u, v, w, e;
        double tol = 0.5 * pow((double)_DBL_RADIX, -DBL_MANT_DIG + 1);
    
        oc = a[0].real();
        od = a[0].imag();
        og = oh = 1.0;
        t = fabs(z.real());
        u = fabs(z.imag());
        for (size_t i = 1; i <= n; i++)
        {
            nc = z.real() * oc - z.imag() * od + a[i].real();
            nd = z.imag() * oc + z.real() * od + a[i].imag();
            v = og + fabs(oc);
            w = oh + fabs(od);
            ng = t * v + u * w + fabs(a[i].real()) + 2.0 * fabs(nc);
            nh = u * v + t * w + fabs(a[i].imag()) + 2.0 * fabs(nd);
            og = ng;
            oh = nh;
            oc = nc;
            od = nd;
        }
        e = abs(complex<double>(ng, nh)) * pow(1 + tol, 5 * n) * tol;
        return e;
    };

多项式归约策略

找到一个根后，我们需要对其进行综合除法以降低多项式次数，为找到下一个根做准备。这时就会出现问题：是使用前向除法还是后向除法？

Wilkinson [7] 表明，为了使归约过程稳定，如果您按模大小递增的顺序找到多项式的根，并且始终先用最小模的根归约多项式，那么您应该选择前向除法；当然，如果您按模大小递减的顺序找到根，则选择后向除法。

为了进行前向除法，我们尝试从最高系数 a_n 开始求解方程

a_nzⁿ+a_n-1z^n-1+⋯+a₁z+a₀=(b_n-1z^n-1+b_n-2z^n-2+⋯+b₁z+b₀)(z-R)

其中 R 是根。

现在解出 bs，得到递推关系

b_n-₁=a_n
b_k=a_k+1+R∙b_k+1, k=n-2,…,0

这个简单的算法适用于实系数和实根的多项式，或者复系数和复根的多项式，只需使用复数算术即可实现相同的递推关系。

// Deflate polynomial and compute new coefficients in coeff
        z = 0;
        for (int j = 0; j < n; j++)
            z = coeff[j] = z * pz.z + coeff[j];
        coeff.resize(n);

K. Madsen 牛顿算法的实现

该求根器的实现遵循 K. Madsen 在 [2] 中首次描述的方法。该方法最初是用 AlgolW 实现的（您还记得那个语言吗？）。下面的实现是修改后的版本，已转换为 C++，并使用了更现代的结构，包括 C++ STL 库。第一步是规划流程。

当然，最有趣的部分是“启动牛顿迭代”部分。Madsen [2] 提供了一个非常快速有效的实现，它不仅在迭代次数不多时就能找到根，而且还能处理牛顿方法通常会遇到的问题。我无意重复 [2] 中已经描述得如此出色的内容，只是想重点介绍他牛顿实现中一些有趣的地方。

首先，我们消除简单根（等于零的根）。
然后我们找到一个合适的起始点来启动牛顿迭代，这还包括基于 P(x) 的可接受值的终止准则，我们将在此停止当前迭代。
启动牛顿迭代
1. 第一步是找到 dz_n=P(z_n)/P’(z_n)，当然，还要决定当 P’(z_n) 为零时应该怎么做。出现这种情况时，大多数时候是因为局部最小值，最好的处理方法是改变方向，乘以因子 dz_n=dz_n(0.6+i0.8)m。这相当于将方向旋转 53 度（对于奇数次幂）并乘以因子 m。当这种情况发生时，m = 5 是一个合理的值。
2. 此外，很容易想到，如果 P’(z_n)~0，您可能会得到一个不合理的步长 dz_n，因此我们引入了一个限制因子，如果 abs(dz_n)>5·abs(dz_n-1)（当前步长大于前一次迭代步长的 5 倍），则减小当前步长。同样，我们通过 dz_n=dz_n(0.6+i0.8)5(abs(dz_n-1)/abs(dz_n)) 来改变方向。
3. 这两个修改（a 和 b）使他的方法非常有弹性，并确保它总是收敛到根。
4. 下一个问题是如何处理重数 > 1 的问题，这会将 2 阶收敛速率降低到线性收敛速率。找到合适的 dz_n 并计算新的后，我们检查 P(z_n+1)>P(z_n)：如果是，我们查看修改后的 z_n+1=z_n-0.5dz_n，如果 P(z_n+1)≥P(z_n)，则我们使用原始的 z_n+1 作为下一次迭代的新的起始点。如果不是，那么我们接受 z_n+1 作为更好的选择，并继续查看新修改的 z_n+1=z_n-0.25dz_n。另一方面，如果新的 P(z_n+1)≥P(z_n)，我们使用前一个 z_n+1 作为下一次迭代的新起始点。如果不是，那么我们假设我们正在接近一个新的鞍点，并且方向会改变为 dz_n=dz_n(0.6+i0.8)，我们将 z_n+1=z_n-mdz_n 作为下一次迭代的新起始点。
  另一方面，如果 P(z_n+1)<=P(z_n)：那么我们正在朝着正确的方向前进，然后我们继续在该方向上进行步进，使用 zn+1=z_n-mdz_n，m=2,..,n，只要 P(z_n+1)<=P(z_n)，并为下一次迭代使用最佳的 m。这个过程的好处是，如果存在重数为 m 的根，那么 m 也是步长的最佳选择，这将即使对于多重根也能保持 2 阶收敛速率。
过程 a-d 继续进行，直到达到停止准则，此时根 z_n 被接受并从多项式中归约。使用归约后的多项式开始新的根搜索。

在 [2] 中，他们将迭代分为两个阶段：阶段 1 和阶段 2。在阶段 1 中，我们试图进入收敛圆，在该圆内我们可以确定牛顿方法将收敛到根。当我们进入该圆时，我们自动切换到阶段 2。在阶段 2 中，我们跳过步骤 d)，仅使用未修改的牛顿步长直到满足停止准则。如果超出收敛圆，我们切换回阶段 1 并继续迭代。

在 [2] 中，他们使用以下准则切换到阶段 2，基于 Ostrowski [3] 的定理 7.1，该定理指出，如果 K 是一个圆，其中心为 w-P(w)/P'(w)，半径为 |P(w)/P'(w)|。

那么，如果我们满足以下两个条件，则保证收敛

以 w 为初始值的牛顿迭代将导致在圆 K 内收敛到零。为了简化计算，我们进行两个替换，因为

并且我们用差分近似代替 p''(w)。

现在我们拥有确定何时切换到阶段 2 所需的所有信息。

C++ 源代码

下面的 C++ 代码使用复数系数的多项式来查找多项式根。如果多项式系数为实数，则可以使用相同的算法。详情请参阅 [1]。

/*
 *******************************************************************************
 *
 *                       Copyright (c) 2023
 *                       Henrik Vestermark
 *                       Denmark, USA
 *
 *                       All Rights Reserved
 *
 *   This source file is subject to the terms and conditions of
 *   Henrik Vestermark Software License Agreement which restricts the manner
 *   in which it may be used.
 *
 *******************************************************************************
*/
/*
 *******************************************************************************
 *
 * Module name     :   Newton.cpp
 * Module ID Nbr   :
 * Description     :   Solve n degree polynomial using Newton's (Madsen) method
 * --------------------------------------------------------------------------
 * Change Record   :
 *
 * Version    Author/Date         Description of changes
 * -------  -------------  ----------------------
 * 01.01     HVE/24Sep2023 Initial release
 *
 * End of Change Record
 * --------------------------------------------------------------------------
*/

// define version string
static char _VNEWTON_[] = "@(#)Newton.cpp 01.01 -- Copyright (C) Henrik Vestermark";

#include <algorithm>
#include <vector>
#include <complex>
#include <iostream>
#include <functional>

using namespace std;
constexpr int       MAX_ITER = 50;

// Find all polynomial zeros using a modified Newton method
// 1) Eliminate all simple roots (roots equal to zero)
// 2) Find a suitable starting point
// 3) Find a root using Newton
// 4) Divide the root up in the polynomial reducing its order with one
// 5) Repeat steps 2 to 4 until the polynomial is of the order of two 
//    whereafter the remaining one or two roots are found by the direct formula
// Notice:
//      The coefficients for p(x) is stored in descending order. 
//      coefficients[0] is a(n)x^n, coefficients[1] is a(n-1)x^(n-1),...,  
//      coefficients[n-1] is a(1)x, coefficients[n] is a(0)
//
static vector<complex<double>> PolynomialRoots
       (const vector<complex<double>>& coefficients)
{
    struct eval { complex<double> z{}; complex<double> pz{}; double apz{}; };
    const complex<double> complexzero(0.0);  // Complex zero (0+i0)
    size_t n;       // Size of Polynomial p(x) 
    eval pz;        // P(z)
    eval pzprev;    // P(zprev)
    eval p1z;       // P'(z)
    eval p1zprev;   // P'(zprev)
    complex<double> z;      // Use as temporary variable
    complex<double> dz;     // The current stepsize dz
    int itercnt;    // Hold the number of iterations per root
    vector<complex<double>> roots;  // Holds the roots of the Polynomial
    vector<complex<double>> coeff(coefficients.size()); // Holds the current 
                                                        // coefficients of P(z)

    copy(coefficients.begin(), coefficients.end(), coeff.begin());
    // Step 1 eliminate all simple roots
    for (n = coeff.size() - 1; n > 0 && coeff.back() == complexzero; --n)
        roots.push_back(complexzero);  // Store zero as the root

    // Compute the next starting point based on the polynomial coefficients
    // A root will always be outside the circle from the origin and radius min
    auto startpoint = [&](const vector<complex<double>>& a)
    {
    const size_t n = a.size() - 1;
    double a0 = log(abs(a.back()));
    double min = exp((a0 - log(abs(a.front()))) / static_cast<double>(n));

    for (size_t i = 1; i < n; i++)
        if (a[i] != complexzero)
        {
            double tmp = exp((a0 - log(abs(a[i]))) / static_cast<double>(n - i));
            if (tmp < min)
                min = tmp;
        }

    return min*0.5;
    };

    // Evaluate a polynomial with complex coefficients a[] at a complex point z and
    // return the result
    // This is the Horner's methods
    auto horner = [](const vector<complex<double>>& a, const complex<double> z)
    {
        const size_t n = a.size() - 1;
        complex<double> fval=a.front();
        eval e;

        for (size_t i = 1; i <= n; i++)
            fval = fval * z + a[i];

        e = { z, fval,abs(fval) };
        return e;
    };

    // Calculate an upper bound for the rounding errors performed in a
    // polynomial with complex coefficient a[] at a complex point z.
    // (Grant & Hitchins test)
    auto upperbound = [](const vector<complex<double>>& a, complex<double> z)
    {
        const size_t n = a.size() - 1;
        double nc, oc, nd, od, ng, og, nh, oh, t, u, v, w, e;
        double tol = 0.5 * pow((double)_DBL_RADIX, -DBL_MANT_DIG + 1);
 
        oc = a[0].real();
        od = a[0].imag();
        og = oh = 1.0;
        t = fabs(z.real());
        u = fabs(z.imag());

        for (size_t i = 1; i <= n; i++)
        {
            nc = z.real() * oc - z.imag() * od + a[i].real();
            nd = z.imag() * oc + z.real() * od + a[i].imag();
            v = og + fabs(oc);
            w = oh + fabs(od);
            ng = t * v + u * w + fabs(a[i].real()) + 2.0 * fabs(nc);
            nh = u * v + t * w + fabs(a[i].imag()) + 2.0 * fabs(nd);
            og = ng;
            oh = nh;
            oc = nc;
            od = nd;
        }
        e = abs(complex<double>(ng, nh)) * pow(1 + tol, 5 * n) * tol;
        return e;
    };

    // Do Newton iteration for polynomial order higher than 2
    for (; n > 2; --n)
    {
        const double Max_stepsize = 5.0; // Allow the next step size to be up to 5 times 
                                         // larger than the previous step size
        const complex<double> rotation = complex<double>(0.6, 0.8);  // Rotation amount
        double r;               // Current radius
        double rprev;           // Previous radius
        double eps;             // The iteration termination value
        bool stage1 = true;     // By default it start the iteration in stage1
        int steps = 1;          // Multisteps if > 1
        vector<complex<double>> coeffprime;

        // Calculate coefficients of p'(x)
        for (int i = 0; i < n; i++)
            coeffprime.push_back(coeff[i] * complex<double>(double(n - i), 0.0));

        // Step 2 find a suitable starting point z
        rprev = startpoint(coeff);      // Computed startpoint
        z = coeff[n - 1] == complexzero ? 
                            complex<double>(1.0) : -coeff[n] / coeff[n - 1];
        z *= complex<double>(rprev) / abs(z);

        // Setup the iteration
        // Current P(z)
        pz = horner(coeff, z);
               
        // pzprev which is the previous z or P(0)
        pzprev.z = complex<double>(0);
        pzprev.pz = coeff[n];
        pzprev.apz = abs(pzprev.pz);

        // p1zprev P'(0) is the P'(0)
        p1zprev.z = pzprev.z;
        p1zprev.pz = coeff[n - 1];       // P'(0)
        p1zprev.apz = abs(p1zprev.pz);

        // Set previous dz and calculate the radius of operations.
        dz = pz.z;                 // dz=z-zprev=z since zprev==0
        r = rprev *= Max_stepsize; // Make a reasonable radius of the 
                                   // maximum step size allowed
        // Preliminary eps computed at point P(0) by a crude estimation
        eps = 6 * n * pzprev.apz * pow((double)_DBL_RADIX, -DBL_MANT_DIG);
       
        // Start iteration and stop if z doesnt change or apz <= eps
        // we do z+dz!=z instead of dz!=0. if dz does not change z 
        // then we accept z as a root
        for (itercnt = 0; pz.z + dz != pz.z && pz.apz > eps && 
                                       itercnt < MAX_ITER; itercnt++)
        {
            // Calculate current P'(z)
            p1z = horner(coeffprime, pz.z);
            if (p1z.apz == 0.0)                 // P'(z)==0 then rotate and try again
                dz *= rotation * complex<double>(Max_stepsize);  // Multiply with 5 
                                                // to get away from saddlepoint
            else
            {
                dz = pz.pz / p1z.pz;  // next dz
                // Check the Magnitude of Newton's step
                r = abs(dz);
                if (r > rprev) // Large than 5 times the previous step size
                {   // then rotate and adjust step size to prevent 
                    // wild step size near P'(z) close to zero
                    dz *= rotation * complex<double>(rprev/r);
                    r = abs(dz);
                }
                rprev = r * Max_stepsize;  // Save 5 times the current step size 
                           // for the next iteration check of reasonable step size
                // Calculate if stage1 is true or false. Stage1 is false 
                // if the Newton converge, otherwise true
                z = (p1zprev.pz - p1z.pz) / (pzprev.z - pz.z);
                stage1 = (abs(z) / p1z.apz > p1z.apz / pz.apz / 4.0) || (steps != 1);
            }

            // Step accepted. Save pz in pzprev
            pzprev = pz;
            z = pzprev.z - dz;      // Next z
            pz = horner(coeff, z); //ff = pz.apz;
            steps = 1;

            if (stage1)
            {   // Try multiple steps or shorten steps depending if P(z) 
                // is an improvement or not P(z)<P(zprev)
                bool div2;
                complex<double> zn;
                eval npz;

                zn = pz.z;
                for (div2 = pz.apz > pzprev.apz; steps <= n; ++steps)
                {
                    if (div2 == true)
                    {  // Shorten steps
                        dz *= complex<double>(0.5);
                        zn = pzprev.z - dz;
                    }
                    else
                        zn -= dz;  // try another step in the same direction

                    // Evaluate new try step
                    npz = horner(coeff, zn);
                    if (npz.apz >= pz.apz)
                        break; // Break if no improvement

                    // Improved => accept step and try another round of step
                    pz = npz;
                    if (div2 == true && steps == 2)
                    {   // To many shorten steps => try another direction and break
                        dz *= rotation;
                        z = pzprev.z - dz;
                        pz = horner(coeff, z);
                        break;
                    }
                }
            }
            else
            {   // calculate the upper bound of error using Grant & 
                // Hitchins's test for complex coefficients
                // Now that we are within the convergence circle.
                eps = upperbound(coeff, pz.z);
            }
        }

        // Check if there is a very small residue in the imaginary part by trying
        // to evaluate P(z.real). if that is less than P(z). 
        // We take that z.real() is a better root than z.
        z = complex<double>(pz.z.real());
        pzprev = horner(coeff, z);
        if (pzprev.apz <= pz.apz)
            pz = pzprev;

        // Save the root
        roots.push_back(pz.z);

       // Deflate polynomial and compute new coefficients in coeff
       z = complex<double>(0);
       for (int j = 0; j < n; j++)
           z = coeff[j] = z * pz.z + coeff[j];
       coeff.resize(n);
    }   // End Iteration

    // Solve any remaining linear or quadratic polynomial
    // For Polynomial with complex coefficients a[],
    // The complex solutions are stored in the back of the roots
    auto quadratic = [&](const std::vector<complex<double>>& a)
    {
        const size_t n = a.size() - 1;
        complex<double> v;

        // Notice a[0] is !=0 since all roots=zero has been captured previously
        if (n == 1)
            roots.push_back(-a[1]/a[0]);
        else
        {
            if (a[1] == complexzero)
            {
                v = sqrt(-a[2] / a[0]);
                roots.push_back(v);
                roots.push_back(-v);
            }
            else
            {
                v = sqrt(complex<double>(1.0)-complex<double>
                        (4.0)*a[0]*a[2]/(a[1]*a[1]));
                if (v.real() < 0)
                    v = (complex<double>(-1.0) - v) * a[1] / 
                        (complex<double>(2.0) * a[0]);
                else
                    v = (complex<double>(-1.0) + v) * a[1] / 
                        (complex<double>(2.0) * a[0]);
                roots.push_back(v);
                roots.push_back(a[2] / (a[0] * v));
            }
        }
        return;
    };

    if (n > 0)
        quadratic(coeff);

    return roots;
}

示例 1

以下是上述源代码工作方式的示例。

For the complex Polynomial: +1x^3+(-13-i1)x^2+(44+i12)x+(-32-i32)

Start Newton Iteration for Polynomial=+1x^3+(-13-i1)x^2+(44+i12)x+(-32-i32)
              Stage 1=&gt;Stop Condition. |f(z)|&lt;3.01e-14
              Start    : z[1]=(0.4+i0.2) dz=(4.31e-1+i2.46e-1) |f(z)|=2.6e+1
Iteration: 1
              Newton Step:  z[1]=(1+i0.7) dz=(-5.91e-1-i4.63e-1) |f(z)|=6.3e+0
              Function value decrease=&gt;try multiple steps in that direction
              Try Step:  z[1]=(2+i1) dz=(-5.91e-1-i4.63e-1) |f(z)|=1.1e+1
                      : No improvement=&gt;Discard last try step
Iteration: 2
              Newton Step:  z[2]=(1.0+i1.0) dz=(-1.83e-2-i2.99e-1) |f(z)|=9.0e-1
              In Stage 2=&gt;New Stop Condition: |f(z)|&lt;4.79e-14
Iteration: 3
              Newton Step:  z[2]=(1.0+i1.0) dz=(4.07e-2+i8.94e-3) |f(z)|=1.8e-2
              In Stage 2=&gt;New Stop Condition: |f(z)|&lt;4.65e-14
Iteration: 4
              Newton Step:  z[4]=(1.000+i1.000) dz=(-6.04e-4-i5.04e-4) |f(z)|=6.3e-6
              In Stage 2=&gt;New Stop Condition: |f(z)|&lt;4.65e-14
Iteration: 5
              Newton Step:  z[8]=(1.0000000+i1.0000000) dz=(2.39e-8-i2.81e-7) 
                            |f(z)|=8.2e-13
              In Stage 2=&gt;New Stop Condition: |f(z)|&lt;4.65e-14
Iteration: 6
              Newton Step:  z[14]=(1.0000000000000+i1.0000000000000) 
              dz=(3.32e-14+i1.53e-14) |f(z)|=7.9e-15
              In Stage 2=&gt;New Stop Condition: |f(z)|&lt;4.65e-14 
              Stop Criteria satisfied after 6 Iterations Final Newton  
              z[14]=(1.0000000000000+i1.0000000000000) dz=(3.32e-14+i1.53e-14) 
              |f(z)|=7.9e-15 Alteration=0% Stage 1=17% Stage 2=83%
              Deflate the complex root z=(0.9999999999999999+i0.9999999999999998)

Solve Polynomial=+(1)x^2+(-12-i2.220446049250313e-16)x+(32+i3.552713678800501e-15) 
directly Using the Newton Method, the Solutions are:
X1=(0.9999999999999999+i0.9999999999999998)
X2=(8.000000000000002-i4.440892098500625e-16)
X3=(3.999999999999999+i6.661338147750937e-16)

示例 2

一个具有重根 z=(1+i) 的相同示例。我们看到每一步都是一个双步，与第一个根的重数为 2 相符。

对于复数多项式：+1x^3+(-10-i2)x^2+(16+i18)x+(-i16)

Start Newton Iteration for Polynomial=+1x^3+(-10-i2)x^2+(16+i18)x+(-i16)
              Stage 1=>Stop Condition. |f(z)|<1.07e-14
              Start    : z[1]=(0.2+i0.2) dz=(2.48e-1+i2.21e-1) |f(z)|=9.1e+0
Iteration: 1
              Newton Step:  z[1]=(0.6+i0.6) dz=(-3.76e-1-i3.54e-1) |f(z)|=2.4e+0
              Function value decrease=>try multiple steps in that direction
              Try Step:  z[1]=(1+i0.9) dz=(-3.76e-1-i3.54e-1) |f(z)|=3.7e-2
                      : Improved=>Continue stepping
              Try Step:  z[1]=(1+i1) dz=(-3.76e-1-i3.54e-1) |f(z)|=1.5e+0
                      : No improvement=>Discard last try step
Iteration: 2
              Newton Step:  z[2]=(1.0+i0.96) dz=(3.68e-4-i3.61e-2) |f(z)|=9.2e-3
              Function value decrease=>try multiple steps in that direction
              Try Step:  z[2]=(1.0+i1.0) dz=(3.68e-4-i3.61e-2) |f(z)|=9.6e-7
                      : Improved=>Continue stepping
              Try Step:  z[2]=(1.0+i1.0) dz=(3.68e-4-i3.61e-2) |f(z)|=9.2e-3
                      : No improvement=>Discard last try step
Iteration: 3
              Newton Step:  z[5]=(1.0002+i1.0000) dz=(1.82e-4+i2.89e-5) |f(z)|=2.4e-7
              Function value decrease=>try multiple steps in that direction
              Try Step:  z[5]=(1.0000+i1.0000) dz=(1.82e-4+i2.89e-5) |f(z)|=8.9e-16
                      : Improved=>Continue stepping
              Try Step:  z[5]=(0.99982+i0.99997) dz=(1.82e-4+i2.89e-5) |f(z)|=2.4e-7
                      : No improvement=>Discard last try step
Stop Criteria satisfied after 3 Iterations
Final Newton  z[5]=(1.0000+i1.0000) dz=(1.82e-4+i2.89e-5) |f(z)|=8.9e-16
Alteration=0% Stage 1=100% Stage 2=0%
              Deflate the complex root z=(0.9999999913789768+i0.9999999957681246)
Solve Polynomial=+(1)x^2+(-9.000000008621024-i1.0000000042318753)x+
                 (8.000000068968186+i8.000000033855004) directly
Using the Newton Method, the Solutions are:
X1=(0.9999999913789768+i0.9999999957681246)
X2=(8.000000000000004-i2.6645352478243948e-15)
X3=(1.0000000086210223+i1.0000000042318753)

结论

本文介绍了一种改进的牛顿方法，该方法最初基于 [2]，提高了寻找复系数多项式根的牛顿方法的效率和稳定性。相同的方法可以轻松应用于实系数多项式（请参阅第二部分）。这是第一部分。第二部分处理了我们只找到实系数多项式根的情况。然而，根仍然可能是复数。第三部分展示了实现高阶方法（例如 Halley）所需的调整。第四部分展示了如何轻松地将 Laguerre 等其他方法纳入同一框架。

可以在 Polynomial roots 找到一个基于 Web 的多项式求解器，它演示了其中许多方法。

参考文献

H. Vestermark. A practical implementation of Polynomial root finders. Practical implementation of Polynomial root finders vs 7.docx (www.hvks.com)
Madsen. A root-finding algorithm based on Newton Method, Bit 13 (1973) 71-75。
A. Ostrowski, Solution of equations and systems of equations, Academic Press, 1966。
Wikipedia Horner’s Method: https://en.wikipedia.org/wiki/Horner%27s_method
Adams, D A stopping criterion for polynomial root finding. Communication of the ACM Volume 10/Number 10/ October 1967 Page 655-658
Grant, J. A. & Hitchins, G D. Two algorithms for the solution of polynomial equations to limiting machine precision. The Computer Journal Volume 18 Number 3, pages 258-264
Wilkinson, J H, Algebraic Processes中的舍入误差, Prentice-Hall Inc, Englewood Cliffs, NJ 1963
McNamee, J.M., Numerical Methods for Roots of Polynomials, Part I & II, Elsevier, Kidlington, Oxford 2009
H. Vestermark, “A Modified Newton and higher orders Iteration for multiple roots.” www.hvks.com/Numerical/papers.html
M.A. Jenkins & J.F. Traub, ”A three-stage Algorithm for Real Polynomials using Quadratic iteration”, SIAM J Numerical Analysis, Vol. 7, No.4, December 1970。

历史

2023年10月11日：参考文献标签不正确
2023年10月10日：初始上传版本缺少一些图片