平方根的C语言实现(三) —

　　版权申明：本文为博主窗户(Colin Cai)原创，欢迎转帖。如要转贴，必须注明原文网址

　　http://www.cnblogs.com/Colin-Cai/p/7223254.html 

　　作者：窗户

　　QQ：6679072

　　E-mail：6679072@qq.com

　　了解了浮点数的存储以及手算平方根的原理，我们可以考虑程序实现了。

　　先实现一个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=;i<;i++));do

        #随机产生bits ~，如果是0，代表测试的数就是0

        #如果不是0，则代表要产生的数二进制可以有多少位

        let bits=RANDOM%

        if [ $bits -eq  ]; then

                x=

                y=

        else

                #产生一个bits位的二进制数x

                x=$({

                        #最高位1

                        echo -n

                        #之后每位随机产生

                        for((j=;j<bits;j++));do

                                let x=RANDOM%

                                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

        fi

done

echo OK

　　测试结果表明，我们的C语言还是可以得到正确的结果的。

　　再来回忆下第一节里讲过的浮点数结构，

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

　对于浮点数a*2ⁿ,

1<=a<2，n为整数，

　　如果n是偶数，

　　那么a*2ⁿ的平方根是sqrt(a)*2^n/2，也满足1<=sqrt(a)<2，n/2是整数；

　　如果n为奇数，

　　那么a*2ⁿ的平方根是sqrt(2*a)*2^(n-1)/2，也满足1<=sqrt(2*a)<2，(n-1)/2是整数。

　　所以此处要用a或者2*a来开平方根，

　　回忆一下浮点数的结构，单精度浮点数的精度是23位。

　　表示的是科学计数法a*2ⁿ的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语言实现(三) ——最终程序实现

最新文章

热门文章