記錄幾個開平方演算法

整數開平方演算法:

本演算法只採用移位、加減法、判斷和迴圈實現，因為它不需要浮點運算，也不需要乘除運算，因此可以很方便地運用到各種晶片上去。

我們先來看看10進位下是如何手工計算開方的。
先看下面兩個算式，
x = 10*p + q  (1)
公式(1)左右平方之後得：
x^2 = 100*p^2 + 20pq + q^2 (2)
現在假設我們知道x^2和p，希望求出q來，求出了q也就求出了x^2的開方x了。
我們把公式(2)改寫為如下格式：
q = (x^2 - 100*p^2)/(20*p+q) (3)
這個算式左右都有q，因此無法直接計算出q來，因此手工的開方演算法和手工除法演算法一樣有一步需要猜值。

我們來一個手工計算的例子：計算1234567890的開方

首先我們把這個數兩位兩位一組分開，計算出最高位為3。也就是(3)中的p，最下面一行的334為餘數，也就是公式(3)中的(x^2 - 100*p^2)近似值

       3
    ---------------
    | 12 34 56 78 90
       9
    ---------------
    |  3 34
下面我們要找到一個0-9的數q使它最接近滿足公式(3)。我們先把p乘以20寫在334左邊：

       3  q
    ---------------
    | 12 34 56 78 90
       9
    ---------------
  6q|  3 34
我們看到q為5時(60+q*q)的值最接近334，而且不超過334。於是我們得到：

       3  5
    ---------------
    | 12 34 56 78 90
       9
    ---------------
  65|  3 34
    |  3 25
    ---------------
          9 56
接下來就是重複上面的步驟了，這裡就不再囉嗦了。 

這個手工演算法其實和10進位關係不大，因此我們可以很容易的把它改為二進位，改為二進位之後，公式(3)就變成了：

q = (x^2 - 4*p^2)/(4*p+q) (4)
我們來看一個例子，計算100(二進位1100100)的開方：

      1  0  1  0
    ---------------
    | 1 10 01 00
      1
    ---------------
 100| 0 10 
    | 0 00 
    ---------------
    |   10 01
1001|   10 01
    ---------------
            0 00
這裡每一步不再是把p乘以20了，而是把p乘以4，也就是把p右移兩位，而由於q的值只能為0或者1，所以我們只需要判斷餘數(x^2 - 4*p^2)和(4*p+1)的大小關係，如果餘數大於等於(4*p+q)那麼該上一個1，否則該上一個0。

下面給出完成的C語言程式，其中root表示p，rem表示每步計算之後的餘數，divisor表示(4*p+1)，通過a>>30取a的最高2位，通過a<<=2將計算後的最高2位剔除。其中root的兩次<<1相當於4*p。程式完全是按照手工計算改寫的，應該不難理解。

unsigned short sqrt(unsigned long a){
  unsigned long rem = 0;
  unsigned long root = 0;
  unsigned long divisor = 0;
  for(int i=0; i<16; i++){
    root <<= 1;
    rem = ((rem << 2) + (a >> 30));
    a <<= 2;
    divisor = (root<<1) + 1;
    if(divisor <= rem){
      rem -= divisor;
      root++;
    }
  }
  return (unsigned short)(root);
}

基於二分法的開平方演算法

#include <stdio.h>
//求絕對值
#define abs(x) (x)>0?(x):(-(x))
 
int  main(){
    double val;
    double eps;
    //value及精度
    printf("Input $value $eps->");
    scanf("%lf %lf",&val,&eps);
    if(val < 0 )
    {
        val*=-1;
    }
    double low,high;
    if(val  < 1)
    {
        low = val;
        high = 1;
    }
    else
    {
        low= 1;
        high = val;
    }
    while(1)
    {
        double mid =  (high + low)/2;
        double tmp =  mid * mid;
        if((abs(tmp-val))< eps)
        {
            printf("%f\n",mid);
            break;
        }
        if(tmp > val)
        {
            high = mid;
        }
        else
       {
            low = mid;
        }
    }
}

神秘的0x5f3759df之卡馬克的開平方演算法

float kamake_sqr(float number) {     

    long i;     
    float x, y;     
    const float f = 1.5F;     
    x = number * 0.5F;     
    y = number;     
    i = *(long *) &y;     
    i = 0x5f3759df - (i >> 1);     
    y = *(float *) &i;     
    y = y * (f - (x * y * y));     
    y = y * (f - (x * y * y));     
    return number * y;     
}     
    
main() {     
    printf("sqr(100)=%f", kamake_sqr(100.0));     
    getch();     
}   

輸出：sqr(100)=9.999964

演算法裡面求平方根一般採用的是無限逼近的方法，比如牛頓迭代法，

比如求5的平方根，選一個猜測值比如2，那麼我們可以這麼算
5/2 = 2.5; 2.5+2/2 = 2.25; 5/2.25 = xxx; 2.25+xxx/2 = xxxx ...

這樣反覆迭代下去，結果必定收斂於sqrt(5)

卡馬克牛就牛在選擇了0x5f3759df 這個開始值，使得迭代的時候收斂速度暴漲，對於Quake III所要求的精度10的負三次方，只需要一次迭代就能夠得到結果。

附加一個小故事：

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

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

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