gnu科学库(GSL)三次多项式误差结果

Question

GSL多项式求解器给了我一个不正确的结果。我尝试使用GSL三次多项式求解器求解以下多项式:

x^3-1.96848e20 x^2+9.07605e28 x+9.07605e28 = 0

Wolframalpha 的结果:

1
• 4.61069 e + 8
• 1.96848 e + 20

在<<p> strong>我的程序的结果:

• 2.30531 e + 08年
• 2.30531 e + 08年
• 1.96848 e + 20

我使用的代码:

   #include <iostream>
    #include <gsl/gsl_poly.h>
    using namespace std;
    int main() {
      double result[3]={0,0,0};
      gsl_poly_solve_cubic(-1.96848e20,9.07605e28,9.07605e28, &result[0], &result[1], &result[2]);
      cout << result[0] << endl;
      cout << result[1] << endl;
      cout << result[2] << endl;
      return 0;
    }

怎么了?

事实上，Wolframalpha也给出了奇怪的答案

更新:Wolframalpha都是错误的(更接近正确答案)，我的代码是错误的。

下面是Python (Numpy库)的相同解决方案

In [11]: np.roots(p)
Out[11]: array([  1.96848000e+20,   4.61068948e+08,  -9.99999998e-01])

Answer 1

在这种情况下，GSL中使用的数值方法不够精确。

使用例如Jenkins-Traub算法(代码在这里)