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

Question

使用浮点数，已知二次公式对b^2>>4ac效果不佳，因为它会产生显着性损失，如其解释here。

我被要求找到一种更好的方法来求解二次方程，我知道有 this 算法。还有其他更有效的公式吗？我怎样才能想出更好的公式？我试图对标准方程进行代数运算，但没有任何结果。

Answer 1

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

在这种情况下，您可以使用 online demo, to get these results:

查看其二次正根的输出

Answer 2

这个答案假设这里主要关注的是关于准确性的鲁棒性，而不是关于中间浮点计算中的上溢或下溢的鲁棒性。该问题表明当使用浮点运算直接应用常用数学公式时，对减法抵消问题的认识，以及解决该问题的技术。

另一个需要考虑的问题是术语 b²-4ac 的准确计算。在以下研究报告中对其进行了详细检查：

William Kahan，“关于没有超精确算术的浮点计算成本”，2004 年 11 月 21 日 (online)

最近对 Kahan 笔记的后续工作着眼于计算两种产品差异的更一般性问题 ab-cd：

Claude-Pierre Jeannerod、Nicolas Louvet、Jean-Michel Muller，“进一步分析 Kahan 算法以准确计算 2 x 2 行列式。” 计算数学，卷。 82，第 284 期，2013 年 10 月，第 2245-2264 页 (online)

这利用了融合乘加运算或 FMA，几乎所有现代处理器（包括 x86-64、ARM64 和 GPU）都可以使用它。它作为标准数学函数 fma() 在 C/C++ 中公开。请注意，在没有硬件支持 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。

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

Numerically stable method for solving quadratic equations

floating-point

numerical-methods