哪个浮点比较更准确,为什么？

Question

我正在尝试使用牛顿方法的不同实现来计算平方根.一个重要的决定是什么时候终止算法.

显然,它不会做使用之差的绝对值y*y和x,其中y是的平方根的当前估计x,因为对于大的值x也可能无法代表具有足够精度的平方根.

所以我应该使用相对标准.天真的我会用这样的东西:

static int sqrt_good_enough(float x, float y) {
  return fabsf(y*y - x) / x < EPS;
}

这看起来效果很好.但是最近我开始阅读Kernighan和Plauger的编程风格元素,他们在第1章中给出了相同算法的Fortran程序,其终止标准(用C翻译)将是:

static int sqrt_good_enough(float x, float y) {
  return fabsf(x/y - y) < EPS * y;
}

两者在数学上是等价的,但是有理由选择一种形式而不是另一种吗？

Answer 1

它们仍然不相同; 底部数学在数学上等同于fabsf(y*y - x) / (y*y) < EPS.我与你看到的问题是,如果y*y溢出(可能因为x是FLT_MAX和y是不幸选择),那么可能永远不会发生终止.以下交互使用双精度数.

>>> import math
>>> x = (2.0 - 2.0 ** -52) * 2.0 ** 1023
>>> y = x / math.sqrt(x)
>>> y * y - x
inf
>>> y == 0.5 * (y + x / y)
True

编辑:作为评论(现已删除)指出,在迭代和终止测试之间共享操作也很好.

EDIT2:两者都可能存在欠正常问题x.该专业规范x,以避免两个极端的并发症.

Answer 2

这两者实际上在数学上并不完全等价，除非你为第一个写 fabsf(y*y - x) / (y*y) < EPS 。（抱歉我原来的评论中有错字）

但我认为关键点是使这里的表达式与牛顿迭代中计算 y 的公式相匹配。例如，如果 y 公式为 y = (y + x/y) / 2，则应使用 Kernighan 和 Plauger 的风格。如果是 y = (y*y + x) / (2*y) 您应该使用 (y*y - x) / (y*y) < EPS。

一般来说，终止标准应该是abs(y(n+1) - y(n))足够小（即小于y(n+1) * EPS）。这就是两个表达式应该匹配的原因。如果它们不完全匹配，则由于缩放不同，终止测试可能会判定残差不够小，而 y(n) 的差值小于浮点误差。结果将是无限循环，因为 y(n) 已停止变化并且永远不会满足终止条件。

例如，以下 Matlab 代码与第一个示例完全相同的牛顿求解器，但它会永远运行：

x = 6.800000000000002
yprev = 0
y = 2
while abs(y*y - x) > eps*abs(y*y)
    yprev = y;
    y = 0.5*(y + x/y);
end

它的C/C++版本也有同样的问题。