14

我正在尝试计算平方根的牛顿方法的不同实现。一个重要的决定是何时终止算法。

显然，使用y*y和x的y平方根的当前估计值之间的绝对差是不行的x，因为对于较大的值，x它可能无法以足够的精度表示其平方根。

所以我应该使用相对标准。我天真地会使用这样的东西：

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

这似乎工作得很好。但最近我开始阅读 Kernighan 和 Plauger 的The Elements of Programming Style，他们在第 1 章中给出了相同算法的 Fortran 程序，其终止标准（用 C 语言翻译）是：

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

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

4

2 回答 2

5

它们仍然不等价。底部的在数学上等价于fabsf(y*y - x) / (y*y) < EPS。我看到的问题是，如果溢出y*y（可能是因为不幸x被选中），那么终止可能永远不会发生。以下交互使用双打。FLT_MAXy

>>> 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：两者都可能有 subnormal 的问题x。专业人员正常化以避免两个极端的x并发症。

于 2012-03-02T13:25:19.227 回答
3

这两者在数学上实际上并不完全等价，除非你为第一个写 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++ 版本也有同样的问题。

于 2012-03-02T14:13:46.310 回答