平方根的C语言实现(一) —— 浮点数的存储

原文链接

 

 曾经做一个硬件成本极度控制的项目,因为硬件成本极低,并且还需要实现较高的精度测量,过程中也自己用C语言实现了正弦、余弦、反正切、平方根等函数。

  以下,无论是在我的实际项目中还是本地的计算机系统,int都是4个字节且机器为小端,除非特别提及,否则都如此默认。按理float的存储没有大小端之分,可是的确在powerpc大端上浮点数的存储也一样是和X86/ARM这样的小端机相反。不过因为正好因大小端而决定浮点数的存储顺序,那么本系列贴子里所有的C语言程序至少在powerpc大端上也是效果相同的。

  尽管在这个项目中我非常想用double来存储小数,但因为这需要翻一倍的存储,从而只好作罢,为了那可怜的存储,我一度甚至想考虑实现3字节的浮点数来,但大致估算了误差(至于如何估算一个公式计算的误差,需要先利用浮点数的结构求自变量的误差,然后要用到数值分析求公式误差,以后有机会开贴单说),实在不靠谱,于是还是用float单精度4字节来存储浮点数。此项目后面围绕着精度、运算时间,从而调整了好几次数据产生、乃至算法的原理,不过这都不是这个系列里要讲的。本系列只讲单精度4字节浮点数的平方根实现,一共分为三节:

  第一节讲浮点数的存储;

  第二节讲手算平方根的原理;

  第三节讲C语言最终实现。

  我们先看浮点数是如何表示实数的,IEEE 754定义了浮点数的结构:

  在了解浮点数的存储之前,我们了解一下科学计数法。

  我们平常用的进位制为十进制,所有不为0的实数都可以表示为s*a*10n,其中:s取1或-1,称为符号部分;a满足1≤a<10,称为小数部分;n为整数,称为指数部分。

  我们的计算机计数一般使用二进制,其道理不用我多说,浮点数也一样用的二进制,用的是二进制下的科学计数法。仿照之前十进制下的科学计数法,即可得二进制下的科学计数。所有不为0的实数,都可以表示为s*a*10n,其中:s取1或-1,称为符号部分;a满足1≤a<2,称为小数部分;n为整数,称为指数部分。32个bit中,最高位1个bits表示符号位s,紧接着的8个bits表示指数位,最后的23个bits表示a。

 

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

 

  用大写表示,代表二进制,与科学计数法的s/n/a关系如下:

  若S为0,则s取1;若S为1,则s取-1。

  n = N(2) - 127,这里N(2)是N的二进制值,

  在符号不会混乱的时候,下面就用N来代替N(2)

  a = 1 + A(2)*2-23 ,这里是A的二进制值,

  在符号不会混乱的时候,下面就用A来代替A(2)

  浮点数和定点数一样,也是离散的,4字节浮点数有32个bits,所以最多只能表示232个不同的实数,是对实数的一种近似,但却有很大的范围,可以满足我们很多的需求了。

  写一个C语言程序来验证这点:

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
#include <stdio.h>
#include <string.h>
#include <inttypes.h>
int main(int argc, char **argv)
{
        union {
                float f;
                uint32_t u;
        } n;
        int i;
 
        scanf("%f", &n.f);
        for(i=31;i>=0;i--) {
                if(n.u&(1u<<i))
                        printf("1");
                else
                        printf("0");
                if(i==31 || i==23)
                        printf(" ");
        }
        printf("\n");
        printf("S:%u\n", (n.u&(1u<<31))>>31);
        printf("N:%u\n", (n.u&(0xff<<23))>>23);
        printf("A:%u\n", n.u&0x7fffff);
        return 0;
}

  随便找几个数验证一下

复制代码
$ echo 1 | ./a.out
0 01111111 00000000000000000000000
S:0
N:127
A:0
$ echo -1 | ./a.out
1 01111111 00000000000000000000000
S:1
N:127
A:0
$ echo 2 | ./a.out
0 10000000 00000000000000000000000
S:0
N:128
A:0
$ echo 3 | ./a.out
0 10000000 10000000000000000000000
S:0
N:128
A:4194304
$ echo 3.5 | ./a.out
0 10000000 11000000000000000000000
S:0
N:128
A:6291456
$ echo 3.75 | ./a.out
0 10000000 11100000000000000000000
S:0
N:128
A:7340032
$ echo 0.75 | ./a.out
0 01111110 10000000000000000000000
S:0
N:126
A:4194304
$ echo 0.875 | ./a.out
0 01111110 11000000000000000000000
S:0
N:126
A:6291456
复制代码

上面的数都是满足的。

可是我们再回头看看,科学计数法其实也是有缺陷的,0其实是无法用科学计数法表示的,可是0使用的场合非常多,所以浮点数必须支持。于是单精度浮点数的2^32种表示中,不是每一种都是科学计数法。

IEEE754规定,单精度浮点数还支持非规格化的数,也就是不是科学计数法的数。

当指数位N为0,也就是N的8个bits全是0的时候,符号位依然是符号位,

表示的数值是s* A*2-149

之所以后面的指数是149,是因为规格化的数所能表示的最小正数为2-126,

而N为0的时候所表示的最大数为(223-1)*2-149,

两者十分接近,

$ echo 'scale=60;2^(-126);(2^23-1)*2^(-149);' | bc
.000000000000000000000000000000000000011754943508222875079687
.000000000000000000000000000000000000011754942106924410159919

编个C语言程序验证一下

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
#include <stdio.h>
#include <string.h>
#include <inttypes.h>
int main(int argc, char **argv)
{
        union {
                float f;
                uint32_t u;
        } n;
        int i;
 
        scanf("%" PRIx32, &n.u);
        for(i=31;i>=0;i--) {
                if(n.u&(1u<<i))
                        printf("1");
                else
                        printf("0");
                if(i==31 || i==23)
                        printf(" ");
        }
        printf("\n");
        printf("S:%u\n", (n.u&(1u<<31))>>31);
        printf("N:%u\n", (n.u&(0xff<<23))>>23);
        printf("A:%u\n", n.u&0x7fffff);
        printf("%.60f\n", n.f);
        return 0;
}

  找几个数验证一下

复制代码
$ echo 0x00000001 | ./a.out
0 00000000 00000000000000000000001
S:0
N:0
A:1
0.000000000000000000000000000000000000000000001401298464324817
$ echo 0x007fffff | ./a.out
0 00000000 11111111111111111111111
S:0
N:0
A:8388607
0.000000000000000000000000000000000000011754942106924410754870
$ echo 0x80000001 | ./a.out
1 00000000 00000000000000000000001
S:1
N:0
A:1
-0.000000000000000000000000000000000000000000001401298464324817
$ echo 0x807fffff | ./a.out
1 00000000 11111111111111111111111
S:1
N:0
A:8388607
-0.000000000000000000000000000000000000011754942106924410754870
$ echo 0x00000000 | ./a.out
0 00000000 00000000000000000000000
S:0
N:0
A:0
0.000000000000000000000000000000000000000000000000000000000000
$ echo 0x80000000 | ./a.out
1 00000000 00000000000000000000000
S:1
N:0
A:0
-0.000000000000000000000000000000000000000000000000000000000000
复制代码

可以看到,有0和-0的区别,浮点数的确就有这么神奇。

另外,IEEE754规定,N等于127,也就是这8个bits全为1的时候也是非规格数,分以下三种情况:

S=0,N=127,A=0时,为正无穷大;

S=1,N=127,A=0时,为负无穷大;

N=127,A≠0是,为NAN(not a number)。

同理,我们还是验证一下:

复制代码
$ echo 0x7f800000 | ./a.out
0 11111111 00000000000000000000000
S:0
N:255
A:0
inf
$ echo 0xff800000 | ./a.out
1 11111111 00000000000000000000000
S:1
N:255
A:0
-inf
$ echo 0x7f800001 | ./a.out
0 11111111 00000000000000000000001
S:0
N:255
A:1
nan
$ echo 0xff800001 | ./a.out
1 11111111 00000000000000000000001
S:1
N:255
A:1
nan
复制代码

 

inf和-inf用于两个实数通过运算产生,因其大小上已经超越浮点数最大可程度以表示的实数,只能用无穷大表示,或者浮点数除0。

而nan则是结果已经不是实数范畴了,比如inf参与了运算,再比如,负数开平方根也会产生nan,这是因为浮点数并不是用于直接表示复数,浮点数并非是要直接模拟一个近似的代数闭包。

posted @ 2018-03-21 12:46  通宵敲代码  阅读(667)  评论(0编辑  收藏  举报