求解二次方程的数值稳定方法

Question

使用浮点,据了解,二次公式不会对于b ^ 2 >> 4AC工作得很好,因为它会产生显着的损失,因为它是解释在这里.

我被要求找到一种更好的方法来解决二次方程,我知道有这种算法.还有其他更好的公式吗？我怎样才能想出更好的配方？我试图用代数方式操纵标准方程,没有任何结果.

Answer 1

这个答案假定这里主要关注的是精度方面的鲁棒性,而不是中间浮点计算中溢出或下溢的鲁棒性.这个问题表明当使用浮点算法直接应用常用的数学公式时,减去消除问题的意识,以及解决它的技术.

另一个需要考虑的问题是术语的准确计算b²-4ac.在以下研究报告中详细检查:

William Kahan,"关于没有超精确算术的浮点计算的成本",2004年11月21日(在线)

最近对Kahan笔记的后续工作研究了计算两种产品差异的更一般问题ab-cd:

Claude-Pierre Jeannerod,Nicolas Louvet,Jean-Michel Muller,"进一步分析Kahan算法,精确计算2 x 2行列式." 计算数学,卷.82,第284号,2013年10月,第2245-2264页(在线)

这使用了融合的乘法 - 加法运算或FMA,它几乎适用于所有现代处理器,包括x86-64,ARM64和GPU.它在C/C++中作为标准数学函数公开fma().请注意,在没有FMA硬件支持的平台上,fma()必须使用仿真,这通常很慢,并且发现一些仿真具有严重的功能缺陷.

FMA计算a*b+c使用完整产品(既不圆形也不以任何方式截断)并在最后应用单个舍入.这允许精确计算两个本机精度浮点数的乘积作为两个本地精度浮点数的未评估总和,而不需要在中间计算中使用扩展精度算法:h = a * b并且l = fma (a, b,- h)其中h+l代表产品a*b 完全.这提供了ab-cd如下的有效计算:

/*
  diff_of_products() computes a*b-c*d with a maximum error < 1.5 ulp

  Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller, 
  "Further Analysis of Kahan's Algorithm for the Accurate Computation 
  of 2x2 Determinants". Mathematics of Computation, Vol. 82, No. 284, 
  Oct. 2013, pp. 2245-2264
*/
double diff_of_products (double a, double b, double c, double d)
{
    double w = d * c;
    double e = fma (-d, c, w);
    double f = fma (a, b, -w);
    return f + e;
}

使用此构建块,如果判别式为正,则可以如下高精度计算二次方程的实根:

/* compute the real roots of a quadratic equation: ax² + bx + c = 0, 
   provided the discriminant b²-4ac is positive
*/
void solve_quadratic (double a, double b, double c, double *x0, double *x1)
{
    double q = -0.5 * (b + copysign (sqrt (diff_of_products (b, b, 4.0*a, c)), b));
    *x0 = q / a;
    *x1 = c / q;
}

在中间计算中没有溢出或下溢的测试用例的广泛测试中,在计算解决方案中观察到的最大误差从未超过3 ulps.

Answer 2

用于自动重新排列浮点表达式以减少舍入错误的Herbie工具通常为解决此类错误提供了一个很好的起点。

在这种情况下，您可以使用在线演示查看其二次方根的输出，以获得以下结果：