平方根的C語言實現(一) —— 浮點數的存儲
曾經做一個硬件成本極度控制的項目，因為硬件成本極低，並且還需要實現較高的精度測量，過程中也自己用C語言實現了正弦、余弦、反正切、平方根等函數。
以下，無論是在我的實際項目中還是本地的計算機系統，int都是4個字節且機器為小端，除非特別提及，否則都如此默認。按理float的存儲沒有大小端之分，可是的確在powerpc大端上浮點數的存儲也一樣是和X86/ARM這樣的小端機相反。不過因為正好因大小端而決定浮點數的存儲順序，那麽本系列貼子裏所有的C語言程序至少在powerpc大端上也是效果相同的。
盡管在這個項目中我非常想用double來存儲小數，但因為這需要翻一倍的存儲，從而只好作罷，為了那可憐的存儲，我一度甚至想考慮實現3字節的浮點數來，但大致估算了誤差（至於如何估算一個公式計算的誤差，需要先利用浮點數的結構求自變量的誤差，然後要用到數值分析求公式誤差，以後有機會開貼單說），實在不靠譜，於是還是用float單精度4字節來存儲浮點數。此項目後面圍繞著精度、運算時間，從而調整了好幾次數據產生、乃至算法的原理，不過這都不是這個系列裏要講的。本系列只講單精度4字節浮點數的平方根實現。
我們先看浮點數是如何表示實數的，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語言程序來驗證這點:
#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語言程序驗證一下
#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，這是因為浮點數並不是用於直接表示復數，浮點數並非是要直接模擬一個近似的代數閉包。
平方根的C語言實現(一) —— 浮點數的存儲