平方根的C語言實現(三) ——最終程式實現

標籤：運算   四捨五入   回憶一下   bsp   測試結果   理論   a.out   進位   簡單   

　　瞭解了浮點數的儲存以及手算平方根的原理，我們可以考慮程式實現了。

　　先實現一個64位整數的平方根，根據之前的手算平方根，程式也不是那麼難寫了。

#include <stdint.h>
uint64_t _sqrt_u64(uint64_t a){        int i;        uint64_t res;        uint64_t remain;        //0的平方根是0,特殊處理一下        if(a == 0ull)                return 0ull;        //找到最高位的1，並且產生平方根結果最高位的1        for(i=62;;i-=2)                if(a&(3ull<<i)) {                        res = 1ull;                        remain = ((a&(3ull<<i))>>i) - 1ull;                        i -= 2;                        break;                }        //根據手算平方根的原理，依次產生各位結果        for(;i>=0;i-=2) {                //右移動兩位，並把a接著的兩位併入remain                remain = (remain<<2)|((a&(3ull<<i))>>i);                if(((res<<2)|1ull) <= remain) {                        //產生新一位的1                        remain = remain - ((res<<2)|1ull);                        res = (res<<1)|1ull;                } else {                        //產生新一位的0                        res <<= 1;                }        }        return res;}

　　其實，可以合在一起寫，代碼會短一些，但效率會低那麼一點點，而且編譯器應該不太容易最佳化。

#include <stdint.h>
uint64_t _sqrt_u64(uint64_t a){        int i;        uint64_t res;        uint64_t remain;        res = remain = 0ull;        for(i=62;i>=0;i-=2) {                remain = (remain<<2)|((a&(3ull<<i))>>i);                if(((res<<2)|1ull) <= remain) {                        remain = remain - ((res<<2)|1ull);                        res = (res<<1)|1ull;                } else {                        res <<= 1;                }        }        return res;}

　　不過，我們不需要這個結果。

　　為了驗證其正確性，我們來寫個C語言的main

#include <stdio.h>#include <stdint.h>#include <inttypes.h>uint64_t _sqrt_u64(uint64_t a);int main(){        uint64_t a, b;        scanf("%" PRIu64, &a);        b = _sqrt_u64(a);        printf("%" PRIu64 "\n",b);        return 0;}

　　我們shell程式測試一下，我們當然不可能測試過每一個64bits的數，這個運算量太大，不現實。我們可以用隨機取一部分來測試。

 

#!/bin/bash#編譯gcc -O2 -s sqrt_u64.c main_sqrt_u64.c -o a.out#隨機測試10000個數for((i=0;i<10000;i++));do        #隨機產生bits 0~64，如果是0，代表測試的數就是0        #如果不是0，則代表要產生的數二進位可以有多少位        let bits=RANDOM%65        if [ $bits -eq 0 ]; then                x=0                y=0        else                #產生一個bits位的位元x                x=$({                        #最高位1                        echo -n 1                        #之後每位隨機產生                        for((j=1;j<bits;j++));do                                let x=RANDOM%2                                echo -n $x                        done                })                #用bc將x轉換成十進位                x=$(echo ‘obase=10;ibase=2;‘"$x" | bc)               #用bc計算x的平方根取整，理論上和我們的C語言計算一致                y=$(echo ‘sqrt(‘"$x"‘)‘ | bc)        fi        #z是我們的C語言計算結果        z=$(echo $x | ./a.out)        #比較，如果不一致，就報錯        if [ $y -ne $z ];then                echo $x $y $z error                exit 1        fidoneecho OK

　　測試結果表明，我們的C語言還是可以得到正確的結果的。

　　再來回憶下第一節裡講過的浮點數結構，

　　S(1bits)  |   N(8bits)  |  A(23bits)

 　  對於浮點數a*2n, 

       1<=a<2，n為整數，

　　如果n是偶數，

　　那麼a*2n的平方根是sqrt(a)*2n/2，也滿足1<=sqrt(a)<2，n/2是整數；

　　如果n為奇數，

　　那麼a*2n的平方根是sqrt(2*a)*2(n-1)/2，也滿足1<=sqrt(2*a)<2，(n-1)/2是整數。

　　所以此處要用a或者2*a來開平方根，

　　回憶一下浮點數的結構，單精確度浮點數的精度是23位。

　　表示的是科學計數法a*2n的a減去1的部分，那麼加上整數1可以用二進位24位表示。

　　於是，我們就想，一個二進位48位或47位長的數，平方根是二進位24位。那麼，我們就可以用一個48位或47位的二進位整數的平方根計算結果的小數部分。

　　nan/inf/-inf以及負數的平方根都是nan，

　　0.0的平方根是0.0,

　　-0.0的平方根是-0.0（可能只是某些庫裡是這樣的），

　　以上都可以在計算的時候特殊化一下。

　　規格數（就是用科學計數法表示的浮點數）的平方根也是規格數，

　　S=0，N=0，A>0代表的是A*2-149，也就是(A*2)*2-150，

　　我們稍微計算一下，可以明白，所有的此類數的平方根都在規格數表示的範圍內。

　　於是，有了以下的代碼。

　　

#include <stdint.h>
static uint32_t _sqrt_(uint64_t a){        int i;        uint64_t res;        uint64_t remain;        res = remain = 0ull;        //之前整數平方根被直接最佳化，我們只需要求47位或者48位整數的平方根        for(i=46;i>=0;i-=2) {                remain = (remain<<2)|((a&(3ull<<i))>>i);                if(((res<<2)|1ull) <= remain) {                        remain = remain - ((res<<2)|1ull);                        res = (res<<1)|1ull;                } else {                        res <<= 1;                }        }        return (uint32_t)res;}float mysqrtf(float f){        union {                float f;                uint32_t u;        } n;        uint32_t N,A;        int _N, i;        uint64_t _A;        n.f = f;        if(n.u == 0x80000000 || n.u == 0x00000000) /* 0.0/-0.0 */                return n.f;        N = (n.u&(0xff<<23))>>23;        if(N==0xff||(n.u&0x80000000)) { /* inf/-inf/nan/  f < 0.0*/                n.u = 0x7fc00000; /* nan */                return n.f;        }        if(N!=0x0) { /* 用科學計數法表示的規格數 */                A = (n.u&0x7fffff)|0x800000;                _N = (int)N - 127;                if(N&0x1) {                        _A = (uint64_t)A<<23;                } else {                        _A = (uint64_t)A<<24;                        _N--;                }        } else { //A*2^(-149)這種表示方式的浮點數                //還是需要找最高位                for(i=22;;i--)                        if(n.u&((0x1)<<i))                                break;                //然後需要移位，要區分奇數和偶數                if(i&0x1) {                        _N = i-149;                        _A = (uint64_t)n.u << (46-i);                } else {                        _N = i-150;                        _A = (uint64_t)n.u << (47-i);                }        }        //小數部分        A = _sqrt_(_A);        //指數部分        N = (uint32_t)(_N/2+127);        //得到結果        n.u = (A&0x7fffff)|(N<<23);        return n.f;}

　　同樣，也寫個測試用的程式，對inf/-inf/nan/0.0/-0.0以及負數不測了，這些很簡單。

#include <stdio.h>#include <stdint.h>#include <math.h>#include <time.h>#include <stdlib.h>#include <inttypes.h>int main(int argc, char **argv){        union {                float f;                uint32_t u;        } n;        uint32_t A,N;        float f,f2;        int i;        srand((unsigned)time(NULL));        //隨機10000個資料        for(i=0;i<10000;i++) {                N = rand()%256;                if(N==255)                        N=254;                A = 0x0;                A |= rand()%256;                A |= (rand()%256)>>8;                A |= (rand()%256)>>16;                n.u = (A&0x7fffff)|(N<<23);                f = sqrtf(n.f);                f2 = mysqrtf(n.f);                printf("%.60f %.60f\n",f,f2);        }        return 0;}

　　結果發現，我們的程式和數學庫裡的sqrtf結果有細微差別。

　　於是，我們決定再加個小東西，就是四捨五入。之前我們用的是47位或者48位元開平方，為了四捨五入，我們需要多一位，於是就用49位或者50位元開平方。

　　修改一下mysqrtf，增加兩位拿去開平方，_sqrt_也動一下。

#include <stdint.h>
static uint32_t _sqrt_(uint64_t a){        int i;        uint64_t res;        uint64_t remain;        res = remain = 0ull;        //之前整數平方根被直接最佳化，我們只需要求49位或者50位整數的平方根        for(i=48;i>=0;i-=2) {//這裡之前是46,改成48                remain = (remain<<2)|((a&(3ull<<i))>>i);                if(((res<<2)|1ull) <= remain) {                        remain = remain - ((res<<2)|1ull);                        res = (res<<1)|1ull;                } else {                        res <<= 1;                }        }        return (uint32_t)res;}

float mysqrtf(float f){        union {                float f;                uint32_t u;        } n;        uint32_t N,A;        int _N, i;        uint64_t _A;        n.f = f;        if(n.u == 0x80000000 || n.u == 0x00000000) /* 0.0/-0.0 */                return n.f;        N = (n.u&(0xff<<23))>>23;        if(N==0xff||(n.u&0x80000000)) { /* inf/-inf/nan/  f < 0.0*/                n.u = 0x7fc00000; /* nan */                return n.f;        }        if(N!=0x0) { /* 用科學計數法表示的規格數 */                A = (n.u&0x7fffff)|0x800000;                _N = (int)N - 127;                if(N&0x1) {                        _A = (uint64_t)A<<25;                } else {                        _A = (uint64_t)A<<26;                        _N--;                }        } else { //A*2^(-149)這種表示方式的浮點數                //還是需要找最高位                for(i=22;;i--)                        if(n.u&((0x1)<<i))                                break;                //然後需要移位，要區分奇數和偶數                if(i&0x1) {                        _N = i-149;                        _A = (uint64_t)n.u << (48-i);                } else {                        _N = i-150;                        _A = (uint64_t)n.u << (49-i);                }        }        //小數部分        A = _sqrt_(_A);        //四捨五入        A = (A+(A&0x1))>>1;        //指數部分        N = (uint32_t)(_N/2+127);        //得到結果        n.u = (A&0x7fffff)|(N<<23);        return n.f;}

　　然後再測，準確無誤。於是我們可以完工了。

平方根的C語言實現(三) ——最終程式實現