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

Question

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

floating-pointnumerical-methods

12

使用浮点数，已知当b^2>>4ac时，二次公式的求解效果不佳，因为这将导致精度损失，具体解释请参见此处。

我被要求找到更好的解决二次方程的方法，我知道有这个算法。是否还有其他更好的公式可以使用？我应该如何提出更好的公式？我尝试对标准公式进行代数变换，但没有任何结果。

- John Keeper

3

您可以使用二的幂次方缩放（例如参见W.Y. Sit的《二次规划》http://www.mmrc.iss.ac.cn/ascm/ascm03/sample.pdf），以及鉴别式的高精度计算（W.Kahan，《浮点运算成本论》http://www.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf）。 - gammatester

2

与Math.SE相关：当a非常小或为0时，解二次方程的数值稳定算法 - Ruslan

2个回答

10

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

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

- ztatlock

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- njuffa · Accepted Answer

本回答假设主要关注精度鲁棒性，而不是中间浮点计算中的溢出或下溢的鲁棒性。问题表明意识到当直接使用浮点算术应用常用的数学公式时出现减法抵消的问题，并介绍了解决方法。

另一个需要考虑的问题是计算项b²-4ac的准确性。该问题在以下研究笔记中得到详细探讨:

William Kahan，“On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic”，2004年11月21日(在线) 对Kahan笔记的最新跟进工作研究了计算两个乘积差ab-cd的更一般问题：

Claude-Pierre Jeannerod、Nicolas Louvet、Jean-Michel Muller，“Further analysis of Kahan's algorithm for the accurate computation of 2 x 2 determinants”Mathematics of Computation，Vol. 82，No. 284，2013年10月，pp. 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。