c++ - 求封闭式多项式根的算法

Question

我正在寻找一种稳健算法(或描述算法的论文)，它可以使用封闭形式找到多项式的根(理想情况下直到第 4 个 debree，但任何东西都可以)解决方案。我只对真正的根源感兴趣。

我第一次尝试求解二次方程涉及到这个(我也有类似风格的三次方程/四次方程代码，但现在让我们关注二次方程):

/**
 *  @brief a simple quadratic equation solver
 *
 *  With double-precision floating-point, this reaches 1e-12 worst-case and 1e-15 average
 *  precision of the roots (the value of the function in the roots). The roots can be however
 *  quite far from the true roots, up to 1e-10 worst-case and 1e-18 average absolute difference
 *  for cases when two roots exist. If only a single root exists, the worst-case precision is
 *  1e-13 and average-case precision is 1e-18.
 *
 *  With single-precision floating-point, this reaches 1e-3 worst-case and 1e-7 average
 *  precision of the roots (the value of the function in the roots). The roots can be however
 *  quite far from the true roots, up to 1e-1 worst-case and 1e-10 average absolute difference
 *  for cases when two roots exist. If only a single root exists, the worst-case precision is
 *  1e+2 (!) and average-case precision is 1e-2. Do not use single-precision floating point,
 *  except if pressed by time.
 *
 *  All the precision measurements are scaled by the maximum absolute coefficient value.
 *
 *  @tparam T is data type of the arguments (default double)
 *  @tparam b_sort_roots is root sorting flag (if set, the roots are
 *      given in ascending (not absolute) value; default true)
 *  @tparam n_2nd_order_coeff_log10_thresh is base 10 logarithm of threshold
 *      on the first coefficient (if below threshold, the equation is a linear one; default -6)
 *  @tparam n_zero_discriminant_log10_thresh is base 10 logarithm of threshold
 *      on the discriminant (if below negative threshold, the equation does not
 *      have a real root, if below threshold, the equation has just a single solution; default -6)
 */
template <class T = double, const bool b_sort_roots = true,
    const int n_2nd_order_coeff_log10_thresh = -6,
    const int n_zero_discriminant_log10_thresh = -6>
class CQuadraticEq {
protected:
    T a; /**< @brief the 2nd order coefficient */
    T b; /**< @brief the 1st order coefficient */
    T c; /**< @brief 0th order coefficient */
    T p_real_root[2]; /**< @brief list of the roots (real parts) */
    //T p_im_root[2]; // imaginary part of the roots
    size_t n_real_root_num; /**< @brief number of real roots */

public:
    /**
     *  @brief default constructor; solves for roots of \f$ax^2 + bx + c = 0\f$
     *
     *  This finds roots of the given equation. It tends to find two identical roots instead of one, rather
     *  than missing one of two different roots - the number of roots found is therefore orientational,
     *  as the roots might have the same value.
     *
     *  @param[in] _a is the 2nd order coefficient
     *  @param[in] _b is the 1st order coefficient
     *  @param[in] _c is 0th order coefficient
     */
    CQuadraticEq(T _a, T _b, T _c) // ax2 + bx + c = 0
        :a(_a), b(_b), c(_c)
    {
        T _aa = fabs(_a);
        if(_aa < f_Power_Static(10, n_2nd_order_coeff_log10_thresh)) { // otherwise division by a yields large numbers, this is then more precise
            p_real_root[0] = -_c / _b;
            //p_im_root[0] = 0;
            n_real_root_num = 1;
            return;
        }
        // a simple linear equation

        if(_aa < 1) { // do not divide always, that makes it worse
            _b /= _a;
            _c /= _a;
            _a = 1;

            // could copy the code here and optimize away division by _a (optimizing compiler might do it for us)
        }
        // improve numerical stability if the coeffs are very small

        const double f_thresh = f_Power_Static(10, n_zero_discriminant_log10_thresh);
        double f_disc = _b * _b - 4 * _a * _c;
        if(f_disc < -f_thresh) // only really negative
            n_real_root_num = 0; // only two complex roots
        else if(/*fabs(f_disc) < f_thresh*/f_disc <= f_thresh) { // otherwise gives problems for double root situations
            p_real_root[0] = T(-_b / (2 * _a));
            n_real_root_num = 1;
        } else {
            f_disc = sqrt(f_disc);
            int i = (b_sort_roots)? ((_a > 0)? 0 : 1) : 0; // produce sorted roots, if required
            p_real_root[i] = T((-_b - f_disc) / (2 * _a));
            p_real_root[1 - i] = T((-_b + f_disc) / (2 * _a));
            //p_im_root[0] = 0;
            //p_im_root[1] = 0;
            n_real_root_num = 2;
        }
    }

    /**
     *  @brief gets number of real roots
     *  @return Returns number of real roots (0 to 2).
     */
    size_t n_RealRoot_Num() const
    {
        _ASSERTE(n_real_root_num >= 0);
        return n_real_root_num;
    }

    /**
     *  @brief gets value of a real root
     *  @param[in] n_index is zero-based index of the root
     *  @return Returns value of the specified root.
     */
    T f_RealRoot(size_t n_index) const
    {
        _ASSERTE(n_index < 2 && n_index < n_real_root_num);
        return p_real_root[n_index];
    }

    /**
     *  @brief evaluates the equation for a given argument
     *  @param[in] f_x is value of the argument \f$x\f$
     *  @return Returns value of \f$ax^2 + bx + c\f$.
     */
    T operator ()(T f_x) const
    {
        T f_x2 = f_x * f_x;
        return f_x2 * a + f_x * b + c;
    }
};

代码很烂，我讨厌所有的阈值。但是对于根在 [-100, 100] 区间的随机方程，这还算不错:

root response precision 1e-100: 6315 cases
root response precision 1e-19: 2 cases
root response precision 1e-17: 2 cases
root response precision 1e-16: 6 cases
root response precision 1e-15: 6333 cases
root response precision 1e-14: 3765 cases
root response precision 1e-13: 241 cases
root response precision 1e-12: 3 cases
2-root solution precision 1e-100: 5353 cases
2-root solution precision 1e-19: 656 cases
2-root solution precision 1e-18: 4481 cases
2-root solution precision 1e-17: 2312 cases
2-root solution precision 1e-16: 455 cases
2-root solution precision 1e-15: 68 cases
2-root solution precision 1e-14: 7 cases
2-root solution precision 1e-13: 2 cases
1-root solution precision 1e-100: 3022 cases
1-root solution precision 1e-19: 38 cases
1-root solution precision 1e-18: 197 cases
1-root solution precision 1e-17: 68 cases
1-root solution precision 1e-16: 7 cases
1-root solution precision 1e-15: 1 cases

请注意，此精度与系数的大小有关，系数的大小通常在 10^6 范围内(因此最终精度远非完美，但可能大部分可用)。然而，如果没有阈值，它几乎毫无用处。

我曾尝试使用多精度算法，这通常效果很好，但往往会拒绝许多根，因为多项式的系数不是多精度的，并且某些多项式无法精确表示(如果在一个二次多项式，它主要要么将它分成两个根(我不介意)或者说没有任何根)。如果我想恢复甚至可能稍微不精确的根，我的代码会变得复杂且充满阈值。

到目前为止，我已经尝试使用 CCmath，但要么我无法正确使用它，要么精度非常差。此外，它在 plrt() 中使用迭代(非封闭形式)求解器。

我尝试过使用 GNU 科学库 gsl_poly_solve_quadratic() 但这似乎是一种幼稚的方法，而且在数值上不是很稳定。

天真地使用 std::complex 数字也被证明是一个非常糟糕的主意，因为精度和速度都可能很糟糕(尤其是三次/四次方程，代码中有很多超越数函数)。

将根恢复为复数是唯一的出路吗？这样就不会遗漏任何根，并且用户可以选择根需要的精确度(从而忽略不太精确根中的小虚部)。

Answer 1

这并没有真正回答你的问题，但我认为你可以改进你所得到的，因为当 b^2 >> ac 时你目前有一个“失去重要性”的问题。在这种情况下，您最终会得到一个类似于 (-b + (b + eps))/(2 * a) 的公式，其中 b 的取消可能会丢失许多来自 eps。

处理此问题的正确方法是对一个根使用二次方根的“正常”方程，对另一个根使用鲜为人知的“替代”或“颠倒”方程。您选择哪条路取决于 _b 的符号。

按照以下这几行更改您的代码应该会减少由此产生的错误。

if( _b > 0 ) {
    p_real_root[i] = T((-_b - f_disc) / (2 * _a));
    p_real_root[1 - i] = T((2 * _c) / (-_b - f_disc));
}
else{
    p_real_root[i] = T((2 * _c) / (-_b + f_disc));
    p_real_root[1 - i] = T((-_b + f_disc) / (2 * _a));
}