C/C++手動實現sqrt()

實現sqrt函數功能

1 二分法

/*
 * function:二分法實現sqrt()
 * author:wanglu
 */
float mysqrt_1(float n)
{
    float left = 0, right = n;
    float mid = 0;
    float last;// 保留上一次的結果

    if(n == 1) return 1;// 特判
    if(n < 0) return 0;// 特判

    do{

        if(mid  > n / mid)// 避免溢出
            right = mid;
        else
            left  = mid;

        last = mid;

        mid = left + (right-left) / 2;// 相較(left + right)/2 能避免溢出

    }while(abs(mid-last) > eps);// 相比使用right-left > eps判斷，這樣更加精確

    return mid;
}

• 執行時間

  n	runtime(ns)
  1	90
  2	404
  3333	570
  999999	710

  ​

2 牛頓迭代法

求出根號a的近似值：首先隨便猜一個近似值x，然后不斷令x等於x和a/x的平均數，迭代個六七次后x的值就已經相當精確了。
例如，我想求根號2等於多少。假如我猜測的結果為4，雖然錯的離譜，但你可以看到使用牛頓迭代法后這個值很快就趨近於根號2了：
( 4 + 2/4 ) / 2 = 2.25
( 2.25 + 2/2.25 ) / 2 = 1.56944..
( 1.56944..+ 2/1.56944..) / 2 = 1.42189..
( 1.42189..+ 2/1.42189..) / 2 = 1.41423..
…

這種算法的原理很簡單，我們僅僅是不斷用(x,f(x))的切線來逼近方程x²-a=0的根。根號a實際上就是x²-a=0的一個正實根，這個函數的導數是2x。也就是說，函數上任一點(x,f(x))處的切線斜率是2x。那么，x-f(x)/(2x)就是一個比x更接近的近似值。代入 f(x)=x²-a得到x-(x²-a)/(2x)，也就是(x+a/x)/2。

/*
 * function:牛頓迭代法實現sqrt()
 * author:wanglu
 */
float mysqrt_2(float n)
{
    float init_value = n;// 牛頓法需要選擇一個初始值，這里使等於n
    float x = init_value;// return value
    float last;// 保留上一次的結果

    do{

        last = x;
        x = (x + n/x)/2;

    }while(abs(x - last) > eps);// 比abs(x-n/x)>eps更精確

    return x;
}

• 執行時間

n	runtime(ns)
1	128
2	170
3333	270
999999	300

​

3 神奇的方法

算法的原理其實不復雜,就是牛頓迭代法,用x-f(x)/f’(x)來不斷的逼近f(x)=a的根。

沒錯，一般的求平方根都是這么循環迭代算的但是卡馬克(quake3作者)真正牛B的地方是他選擇了一個神秘的常數0x5f3759df 來計算那個猜測值，就是我們加注釋的那一行，那一行算出的值非常接近1/sqrt(n)，這樣我們只需要2次牛頓迭代就可以達到我們所需要的精度。好吧如果這個還不算NB,接着看:

普渡大學的數學家Chris Lomont看了以后覺得有趣，決定要研究一下卡馬克弄出來的這個猜測值有什么奧秘。Lomont也是個牛人，在精心研究之后從理論上也推導出一個最佳猜測值，和卡馬克的數字非常接近, 0x5f37642f。卡馬克真牛，他是外星人嗎？

傳奇並沒有在這里結束。Lomont計算出結果以后非常滿意，於是拿自己計算出的起始值和卡馬克的神秘數字做比賽，看看誰的數字能夠更快更精確的求得平方根。結果是卡馬克贏了… 誰也不知道卡馬克是怎么找到這個數字的。

最后Lomont怒了，采用暴力方法一個數字一個數字試過來，終於找到一個比卡馬克數字要好上那么一丁點的數字，雖然實際上這兩個數字所產生的結果非常近似，這個暴力得出的數字是0x5f375a86。

/*
 * function:神奇的算法實現sqrt()
 * author:wanglu
 */
float mysqrt_3(float x)
{

    float xhalf = 0.5f*x;
    int i = *(int*)&x; 
    
    if(!x) return 0;
    
    i = 0x5f375a86- (i>>1); // beautiful number
    x = *(float*)&i; 
    x = x*(1.5f-xhalf*x*x); // 牛頓迭代法，提高精度
    x = x*(1.5f-xhalf*x*x); // 牛頓迭代法，提高精度
    x = x*(1.5f-xhalf*x*x); // 牛頓迭代法，提高精度

    return 1/x;
}

• 執行時間

n	runtime(ns)
1	100
2	100
3333	100
999999	100

可以看出這種方法的效率特別高，在實際應用中使用它效果會比其他的好。

​

本文的project

Github倉庫：https://github.com/lrw998/mysqrt

項目使用cmake管理，源文件為mysqrt.cpp

​

使用說明：

（1）cd進入mysqrt目錄

（2）執行cmake .命令生成makefile

（3）執行make編譯工程，在mysqrt/bin目錄下生成可執行文件

（4） cd進入./bin目錄，./mysqrt運行可執行文件

​

引用

https://www.2cto.com/kf/201206/137256.html