平方根的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,这是因为浮点数并不是用于直接表示复数,浮点数并非是要直接模拟一个近似的代数闭包。