平方根算法 C++

Question

我无法弄清楚为什么这个算法在输入的数字超过12位数时进入了一个无限循环。有人能看出来它为什么永远不会结束吗？谢谢。 我刚刚更新了算法，使用了fabs()函数，但仍然进入了一个无限循环。

double squareroot(double x)

{ /* computes the square root of x */

/* make sure x is not negative .. no math crimes allowed! */
assert( x >= 0 );
if (x==0) return 0;

/* the sqrt must be between xhi and xlo */
double xhi = x;
double xlo = 0;
double guess = x/2;

/* We stop when guess*guess-x is very small */

while (abs(guess*guess-x) > 0.00001 )
{
    if (guess*guess > x){
        xhi = guess;
    }

    else {
        xlo = guess;
    }

    guess = (xhi + xlo)/2;
}
return guess;
}

- csciXAV_12

---

abs 是针对整数的函数。你是不是想说 fabs？ - M Oehm

guess*guess 的结果会有12位数字吗？那可真的是很大了。 - Ankur

您的终止准则中的固定_epsilon_值可能也应与输入值相关。否则，例如16.0e-8的计算甚至不会进入第一次迭代。 - M Oehm

@M Oehm 这对于 C 是正确的，但在 C++ 中，<cmath> 也有一个浮点数的重载版本用于 abs。 - TNA

@TNA，是的，我已经承认这一点了。我不是C++程序员，只是通过 _algorithm_来到这里的。我看到的代码是C。 - M Oehm

Answer 1

我认为您应该使用相对误差来终止，而不是绝对误差。

while (abs((guess*guess-x) / guess) > 0.00001)

否则，计算非常长的值的平方根将需要很长时间（它不是无限循环）。

http://en.wikipedia.org/wiki/Approximation_error

干杯！

编辑： 此外，正如下面评论中指出的那样，值得检查是否已经猜过guess，以避免在某些具体情况下陷入无限循环。

- ale64bit

---

当数字很大时，循环可能会变得无限：当 xhi 和 xlo 相邻时，中点将始终是它们之一。当 xhi - xlo 大于所选的 epsilon 时，您将得到一个无限循环。 - M Oehm

@MOehm 不错的观点。这可以进行额外的检查，即当新猜测已经在之前被猜过时应该退出循环（guess == xhi || guess == xlo）。 - leemes

@MOehm 同意。TonyD的答案解决了这个问题，正如leemes在他的评论中指出的那样。（虽然得到一个实际的x来实现这一点会很酷:D） - ale64bit

Answer 2

我建议等待稳定的答案，而不是调整epsilon值：

double squareroot(double x)
{
    if (x < 1) return 1.0 / squareroot(x);  // MSalter's general solution

    double xhi = x;
    double xlo = 0;
    double guess = x/2;

    while (guess * guess != x)
    {
        if (guess * guess > x)
            xhi = guess;
        else
            xlo = guess;

        double new_guess = (xhi + xlo) / 2;
        if (new_guess == guess)
            break; // not getting closer
        guess = new_guess;
    }
    return guess;
}

- Tony Delroy

---

非常好！因为重复使用相同的猜测而导致无限循环的解决方法当然是强制使用不同的猜测。 - M Oehm

这对x=0.25有效吗？因为该算法似乎假定sqrt(x)在范围[0，x]内，而sqrt(0.25)>0.25。第一次猜测是0.125，第二次猜测是0.1875等等。看起来xlo收敛于xhi，并且当xlo=xhi=guess=0.25时循环终止，这显然不是0.5。 - MSalters

@MSalters：不是的...这值得一提；我并不试图修复算法 - 只是在处理问题中的收敛/精度问题。 - Tony Delroy

好的，修复起来很容易。sqrt(1/x) = 1/sqrt(x)，所以只需添加 if (x<1) return 1.0/squareroot(x); 即可。 - MSalters

Answer 3

这不是对你问题的直接回答，而是另一种解决方案。

你可以使用牛顿迭代法求根：

assert(x >= 0);
if (x == 0)
    return 0;

double guess = x;
for (int i=0; i<NUM_OF_ITERATIONS; i++)
    guess -= (guess*guess-x)/(2*guess);
return guess;

24次迭代应该可以得到足够好的近似值，但你也可以检查绝对差异。

- barak manos

Answer 4

我认为当数字足够大时，您无法使用绝对epsilon值，因为它不适合精度。

尝试改用相对比较。请考虑以下函数以检查2个双精度浮点数是否相等：

bool are_eql(double n1, double n2, double abs_e, double rel_e)
{
  // also add check that n1 and n2 are not NaN
  double diff = abs(n1 - n2);
  if (diff < abs_e)
    return true;
  double largest = max(abs(n1), abs(n2));
  if (diff < largest * rel_e)
    return true;
  return false;
}

- denis-bu

---

第一行显然应该是 abs(n1-n2)。我猜 largest 应该是 abs(max(n1,n2)) - 如果两个数都是负数怎么办？在你的代码中，largest 可能是负数，在这种情况下 diff>largest（因为 +0.01 > -100.0）。 - MSalters

Answer 5

http://floating-point-gui.de/errors/comparison/

嘿，看起来你不应该那样使用abs()。有些情况下它应该停止，但由于精度有限，它不会停止。

相反，请使用fabs()

http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

- Álvaro Gómez

---

谢谢，我刚刚更新了使用fabs()函数的代码，但仍然出现无限循环。 - csciXAV_12

似乎有一些C++程序员不知道std::abs函数的重载用法。 - n. m.

@n.m.：说得好。我是通过“算法”标签来到这里的，没有看到“C++”标签。只是看了代码，发现它是用C写的。没有流、命名空间和<math>。 - M Oehm

@csciXAV_12 std::abs有一个重载版本（可能还有头文件<cmath>中的abs）。请参阅http://en.cppreference.com/w/cpp/numeric/math/fabs。发布一个MCVE很重要，如果可能的话，请明确使用`std::abs`。 - juanchopanza