1
2
3
4
5
6
7
8
9
版权申明:本文为博主窗户(Colin Cai)原创,欢迎转帖。如要转贴,必须注明原文网址
 
http://www.cnblogs.com/Colin-Cai/p/7223254.html
 
作者:窗户
 
QQ:6679072
 
E-mail:6679072@qq.com

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

  先实现一个64位整数的平方根,根据之前的手算平方根,程序也不是那么难写了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#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;
}

  其实,可以合在一起写,代码会短一些,但效率会低那么一点点,而且编译器应该不太容易优化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#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

1
2
3
4
5
6
7
8
9
10
11
12
13
#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
        fi
done
echo 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

  我们稍微计算一下,可以明白,所有的此类数的平方根都在规格数表示的范围内。

  于是,有了以下的代码。

  

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#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以及负数不测了,这些很简单。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#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_也动一下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#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;
}<br>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;
}

  然后再测,准确无误。于是我们可以完工了。

posted on   窗户  阅读(7479)  评论(0编辑  收藏  举报
编辑推荐:
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
· 没有源码,如何修改代码逻辑?
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
阅读排行:
· 在鹅厂做java开发是什么体验
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
< 2025年2月 >
26 27 28 29 30 31 1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 1
2 3 4 5 6 7 8

点击右上角即可分享
微信分享提示