Like World

Like's blog

求平方根算法

1. 牛顿法

先猜一个初值，然后依次迭代，直到结果足够精确。比如求number的平方根：

  • X0 = guess
  • X1 = (X0 + number / X0) / 2
  • X2 = (X1 + number / X1) / 2
  • ……
  • Xn = (Xn-1 + number / Xn-1) / 2

结束条件可以是 abs(Xn * Xn – number) 小于某个数，或者两次迭代结果 abs(Xn – Xn-1) 小于某个数.

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

#include <iostream>
#include <cmath>
#include <limits>

using namespace std;

double sqrt_newton(double x)
{
    if (x <= 0)
        return x;

    double val = x, prev;

    do {
        prev = val;
        val = (prev + x / prev) / 2;
    } while(abs(val - prev) > numeric_limits<double>::epsilon());

    return val;
}

注意：上面这个代码不能使用第一个条件判断。因为此时求sqrt_newton(3)时，abs(val * val – x)最后一直是4.44089e-16，而numeric_limits::epsilon()是2.22045e-16。这样就导致了死循环。这应该是由于浮点数的精度问题导致的。 如果要使用第一个条件，可以用 abs(val * val – x) > 0.00000001。

2. 二分法

二分法就是每次取中间的值，然后比较。直到结果足够精确。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

double sqrt_bi(double x)
{
    if (x <= 0)
        return x;

    double low = 0, high = x;
    double mid = (low + high) / 2, prev;

    do {
        if (mid * mid > x)
            high = mid;
        else
            low = mid;

        prev = mid;
        mid = low + (high - low) / 2;
    } while(abs(mid - prev) > numeric_limits<double>::epsilon());

    return mid;
}

3. 不知所云法

QUAKE-III里用来计算1/sqrt(x)的算法，速度快。 这里的float不能改成double，改成double就不对了。

还有个问题是，下面这个算法求出的结果精度会有点问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

float sqrt_magic(float x)
{
    long i;
    float x2, y;
    const float threehalfs = 1.5;

    x2 = x * 0.5;
    y = x;
    i = * (long *) &y;
    i = 0x5f3759df - (i >> 1);
    y = * (float *) &i;
    y = y * (threehalfs - (x2 * y * y));
    y = y * (threehalfs - (x2 * y * y));

    return 1 / y;
}

补充知识：C++中的numeric_limits

numeric_limits定义于limits头文件中，主要用于取代C语言中limits.h中的那些宏。具体见它的参考页面

比如，

  • INT_MAX和INT_MIN使用numeric_limits<int>::max()和numeric_limits<int>::min()代替。

  • DBL_EPSILON使用numeric_limits<double>::epsilon()代替

References：

[1] 一个Sqrt函数引发的血案

[2] 维基百科-牛顿法