floating-point - 为什么 (1-x)(1+x) 比 (1-x^2) 更受欢迎？

Question

我正在查看 arcsin 的运行时库实现这是通过计算实现的:

ArcTan(X, Sqrt(1 - X*X))

但是，计算 1 - X*X的代码实际评估 (1-X)*(1+X) .选择后者有充分的理由吗？我怀疑后者会降低 X 的舍入误差接近于零，但我无法解释为什么会这样。

Answer 1

ArcTan(X, Sqrt(1-X*X))的导数相对于 X 是 1/Sqrt(1-X*X) .这趋于无穷大为 |X|变为 1。因此，当 X 接近 1 或 -1 时，评估中的任何错误都会对结果产生巨大影响。因此，在这些情况下，评估最小化错误是至关重要的。
当 X 接近 1 时，1-X 的评估没有错误(在 IEEE 754 或任何好的浮点系统中，因为结果的标度使得其最低有效位至少与 1 或 X 中的最低有效位一样低，因此精确的数学结果没有可用有效位之外的位)。自 1-X是准确的，请考虑 1+X 中错误的影响通过考虑 ArcTan(X, Sqrt((1-X)*(1+X+e)) 的导数关于 e，其中 e 是 1+X 中引入的误差手术。当 X 接近 1 且 e 很小时，导数约为 -1/10。 (用 Maple 求导并用 1 代替 x 产生 -1/(sqrt(4+2e)*(5+2e) 。然后用 0 代替 e 产生 -1/10。)因此，1+X 中的误差不重要。
因此，将表达式计算为 ArcTan(X, Sqrt((1-X)*(1+X))是一个很好的评价方式。
对于接近 -1 的 X，情况是对称的。 ( 1+X 没有错误，1-X 不是关键。)
相反，如果我们考虑 X*X 中的错误，ArcTan(X, Sqrt(1-X*X+e)) 的导数关于e是，当X接近1时，大约是-1/(2sqrt(e)(1+e))，所以当e小时它就大。所以在评估 X*X 时有一个小错误当 X 接近 1 时，会导致结果出现较大误差。

Ask Pascal Cuoq 指出，在评估函数 f(x) 时，我们通常对最小化最终结果中的相对误差感兴趣。而且，正如我所指出的，计算过程中出现的错误通常是由于浮点舍入导致的中间结果中的相对错误。我能够在上面忽略这一点，因为我正在考虑 X 接近 1 时的函数，因此所考虑的中间值(1+X 和 X*X)和最终值的幅度都接近 1，因此将这些值相除那些幅度不会有任何显着变化。
但是，为了完整起见，我更仔细地检查了情况。在Maple，我写了g := arctan(x, sqrt((1-x*x*(1+e0))*(1+e1))*(1+e2)) ，因此在 x*x 的计算中允许相对误差 e0、e1 和 e2 , 1-x*x ，以及 sqrt分别和我写了h:= arctan(x, sqrt((1-x)*(1+x)*(1+e0))*(1+e2))对于替代方案。请注意，在这种情况下，e0 结合了 1-x 中的三个错误。 , 1+x ，以及它们的乘法；完整的错误项可能是 (1+ea)*(1+eb)*(1+ec) ，但这实际上是 1+e0 e0 的可能范围更大。
然后我检查了这些函数关于(一次一个)e0、e1 和 e2 除以 abs(f(x)) 的导数，其中 f是理想的功能，arctan(x, sqrt(1-x*x)) .例如，在 Maple 中，我检查了 diff(g, e0) / abs(f(x)) .我没有对这些进行全面的分析评估；我检查了一些接近 0 和接近 1 的 x 值以及 e0、e1 和 e2 在它们的极限之一 -2-54 处的值。
对于接近 0 的 x，这些值的大小都约为 1 或更小。也就是说，计算中的任何相对误差都会导致结果中的相对误差类似或更少。
对于接近 1 的 x，e1 和 e2 的导数的值很小，大约为 10-8 或更小。然而，这两种方法的 e0 导数的值有很大的不同。对于 1-x*x方法，该值约为 2•107(使用 x = 1-2-53)。对于 (1-x)*(1+x)方法，该值约为 5•10-9。
综上所述，这两种方法在 x = 0 附近差别不大，但 (1-x)*(1+x)方法在 x = 1 附近明显更好。