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

原文链接

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

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

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

　　第一节讲浮点数的存储；

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

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

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

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

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

　　我们的计算机计数一般使用二进制，其道理不用我多说，浮点数也一样用的二进制，用的是二进制下的科学计数法。仿照之前十进制下的科学计数法，即可得二进制下的科学计数。所有不为0的实数，都可以表示为s*a*10ⁿ，其中：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，所以最多只能表示2³²个不同的实数，是对实数的一种近似，但却有很大的范围，可以满足我们很多的需求了。

　　写一个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(" ");         printf("S:%u ", (n.u&(1u<<31))>>31);         printf("N:%u ", (n.u&(0xff<<23))>>23);         printf("A:%u ", 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的时候所表示的最大数为(2²³-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(" ");         printf("S:%u ", (n.u&(1u<<31))>>31);         printf("N:%u ", (n.u&(0xff<<23))>>23);         printf("A:%u ", n.u&0x7fffff);         printf("%.60f ", 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，这是因为浮点数并不是用于直接表示复数，浮点数并非是要直接模拟一个近似的代数闭包。