平方根的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(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，這是因為浮點數並不是用於直接表示複數，浮點數並非是要直接模擬一個近似的代數閉包。