如何计算符号（a² -b * c） * sqrt（abs（a² -b * c））准确地使用浮点算术？

Question

问题 我想知道是否有一种计算数字准确的方法

sign(a² - b * c) * sqrt(abs(a² - b * c))

浮点算术，因为它遭受了（从大多数到最小问题）：

cancellation n时cancellation n时，这同时使

a² ~= b * c

-和
• sign

  因子不稳定

  由于方形/乘积在两个术语之间放大任何差异流flow/下水流问题（即使您需要很大的数字才能发生，所以我认为这不是一个特别大的问题）
如果可能的话，我想坚持使用32位和64位浮点（例如较高的精度数据类型）（例如小数）或任意精确库（例如Python
• sqrt

  ）等公共浮点数据类型。
  background

  当研究

  a²
问题
• 时，我找到了代码的一部分，看起来在数字上不稳定（
link

），甚至在评论中均表明：

b * c

该代码基本上是针对使用

计算

mpmath

`rnorm->a`

dampsq->b`

`Xxnorm->c`

• 到目前为止我尝试了什么

  这个问题看起来与存在补偿快速算法的准确计算相似，例如，如所提出

  Borgesc.f.，
快速补偿的算法，用于互惠平方根，相互斜侧和givens旋转
• ，

  arxiv：2103.08694

• 计算利用融合的多元ADD操作（FMA）并补偿浮点错误。

  ，无论如何，我对浮点数数学数学的这种深层理论并不特别熟悉，因此我无法将算法转换为

  scipy
的计算。

看来，问题的问题与

减法而不是加法

绝对术语

符号前firactor

        # Distinguish between
        #    r1norm = ||b - Ax|| and
        #    r2norm = rnorm in current code
        #           = sqrt(r1norm^2 + damp^2*||x - x0||^2).
        #    Estimate r1norm from
        #    r1norm = sqrt(r2norm^2 - damp^2*||x - x0||^2).
        # Although there is cancellation, it might be accurate enough.
        if damp > 0:
            r1sq = rnorm**2 - dampsq * xxnorm
            r1norm = sqrt(abs(r1sq))
            if r1sq < 0:
                r1norm = -r1norm

是一个数字问题。

没有诉诸更广泛的类型，一个改进时，当：确定

sign(a² - b * c) * sqrt(abs(a² - b * c))

然后形成产品

hypot(a, b) = sqrt(a² + b²)

。

代码也可以使用其他技巧，如果
• range
也很重要，但听起来像是OP担心精度。 我已经不时地将其与二次方程式形式一起使用。