今天看到一個浮點運算的問題, 覺得很有代表性, 就記錄一下. 有同事問下面這段代碼運行結果是什么, 還很貼心的給了4個選項:
A. 30500 (0.000061 * 500 * 1000000)
B. 與A同量級的近似值
C. 具體值不知道, 但肯定是小於A量級某個浮點數
D. 1024.0 (精確值)
#include <stdlib.h>
#include <stdio.h>
#define LOOP 1000000
#define NUM 500
#define DATA 0.000061
int main() {
float *pf = malloc(NUM * sizeof(float));
for (int i = 0; i < NUM; ++i)
pf[i] = DATA;
float res = 0.0;
for (int j = 0; j < LOOP; ++j)
for (int i = 0; i < NUM; ++i)
res += pf[i];
printf("%f\n", res);
return 0;
}
問了周圍幾個人都回答B或者C, 其中答B居多. 能夠答C的已經考慮比較周全了, 但是還是對IEEE754的格式沒有明確的理解.
先說下答案是D, 實際上在循環次數足夠大的情況下其結果僅與DATA值相關. 為什么? 這與IEEE754的編碼與計算方式有關.
浮點數的由來與IEEE754格式
早期的處理器只支持整型運算單元, 每個整型數據一對一唯一的編碼成二進制數據(i.e. 5 = 0b101, 11 = 0b1011), 制約數值表達范圍的因素僅有一個: 數據編碼長度.
當然編碼長度與硬件成本息息相關, 不同的變量有不同的取值范圍, 所以程序員設計了char, short, int, long等不同位寬的數據類型來支持不同范圍的數據運算.
但是當編碼擴展到小數時問題就來了: 編碼的最低位等於整數1(即整型數據的表示方式的最小間距是1), 那該如何表示一個小於1的間距(值)呢?
基礎的想法是像自然數學一樣擴展, 以某個位為起始, 其低位的位用來表示小數部分, 離基礎位越近的權重越大. i.e. 若取低8位作為小數部分則有
1 = 0x100
1000 = 0x3e800
0.5 = 1 / 2 = 0x100 >> 1 = 0x80
0.125 = 0x20
1000.125 = 0x3e820
這種定點數的編碼方式有幾個顯而易見的問題:
- 表達的精度有限, 誤差為2的-n次(n為小數編碼位數). 如果想要擴展精度只有小數部分的編碼位數.
- 由於小數部分占據的編碼低位, 不同精度的定點數不能直接通過加法器相加(高精度的編碼的小數部分的MSB是低精度的編碼的整數部分的LSB).
- 當編碼的數值較大以至於小數部分可以忽略時(i.e. 上文中的1000.125)會造成編碼空間的浪費(8bit編碼僅表達數值的萬分之一).
為解決以上問題程序員將科學計數法引入編碼中, 現在一個數值的編碼由兩部分相乘組成: 以2為底的指數(被稱為exponent)與一個小於1的分數(被稱為fraction).
在這種編碼方式下, 編碼的一部分, 比如高8位作為指數編碼空間, 剩余的空間用來編碼分數.
4 = 8 * 0.5 = pow(2, 3) * 0.5 = 0x03 << 24 | 0x800000 = 0x03800000
12 = 16 * 0.75 = 0x04c00000
以4為例首先將其轉為科學計數法(pow(2, 3) * 0.5), 然后分別對指數部分與分數部分編碼, 最后拼接在一起得到最終值.
這種編碼方式下分數的每位的權重(對應實際的數值)是非固定的(由指數位決定), 因此也被稱為浮點數.
IEEE754是一種實現浮點數及其運算的標准, 它被廣泛使用在許多現代處理器架構中, 其文檔可以參見這里.
IEEE754的編碼方式與上文稍稍有些不同, 其編碼公式為(pow(2, exponent - bias) * (1 + fraction)), 這里以常見的單精度(float)與雙精度(double)為例:
| type | length | exponent | fraction |
|---|---|---|---|
| float | 32 | 8 | 24 |
| double | 64 | 11 | 53 |
如上表所示float類型的指數位8位, 即指數范圍0 - 255. 為支持小於1的數值編碼還需要減去一個指數偏移值(bias), 其大小為(pow(2, exponent - 1) - 1)(取一半再減1). 因此float實際所能表達的指數范圍是-127 - 128.
float類型的分數位24位, 其中1位符號位(占據編碼最高位). 因此float類型的精度是2的-23次. 注意IEEE754規定了fraction僅包含小數部分, 即運算結果還需要加1.
從網上拖的圖為例:
0.15625 = 1.25 / 8 = pow(2, -3) * 1.25 = pow(2, 124 - 127) * (1 + 0.25) = 124 << 23 | 0x200000 = 0x3e200000
首先將0.15625轉換為科學計數法, 得到指數部分為(1 / 8), 分數部分(1 + 0.25). 由於指數需要減去偏移值, 所以被編碼的指數值為(-3 + 127) = 124, 而分數部分被編碼值為0.25, 最終拼接值為0x3e200000.
寫個case證明一下:
[23:33:52] hansy@hansy:~/testcase$ cat 1.c
#include <stdlib.h>
#include <stdio.h>
int main() {
float a = 0.15625;
int *pa = &a;
printf("%x\n", *pa);
a = 0.25;
printf("%x\n", *pa);
return 0;
}
[23:34:00] hansy@hansy:~/testcase$ gcc 1.c -w && ./a.out
3e200000
3e800000
浮點數的算術運算
IEEE754要求浮點數在做算術運算時需要首先對齊到同一指數位, 然后再做運算, 對齊后做運算的過程會丟失精度. 讓我們來看個例子:
[00:47:18] hansy@hansy:~/testcase$ cat 1.c
#include <stdlib.h>
#include <stdio.h>
int main() {
float a = 1.00001;
float b = 100.0;
int *pa = &a;
int *pb = &b;
printf("%x\n", *pa);
printf("%x\n", *pb);
a += b;
printf("%x\n", *pa);
printf("%f\n", a);
return 0;
}
[00:47:21] hansy@hansy:~/testcase$ gcc 1.c -w && ./a.out
3f800054
42c80000
42ca0001
101.000008
1.00001與100.0相加, 結果不是101.00001, 也不是101.00000, 而是101.000008, 這是為什么呢?
1.00001對應的指數位為0(指數編碼127), 分數位為0x54(84 / (1 << 23)). 100.0對應的指數位為6(指數編碼133), 分數位為0x480000(1.5625 = 1 + 0.5625 = 1 + (0.5 + 0.0625)).
由於兩者量級不同, 1.00001首先會對階到同一指數, 1.00001指數位增加6, 因此對應的分數部分右移6位(指數權重擴大需要對應分數位降低權重), 相當與變成了0x42800000 | 0x1.
然后分數部分兩者相加0x480000 + 0x1 = 0x480001, 注意前面計算的1.00001的分數部分並不包含那個1, 實際上移位后1也需要轉換到對應的分數編碼0x20000(指數位為6, 需要右移6次, 即1 << (23 - 6)), 所以最終分數編碼為0x4a0001. 轉換為float數值為101.00008
從這個例子可以看出, 浮點運算在做對階操作過程中會產生精度的丟失. 當參與運算的兩個數指數相差很大, 大到超過分數編碼長度時, 對階后分數部分為0.
回到開頭
為什么最終答案是1024? 因為float分數編碼23+1位(23位分數編碼加1對應的整型位), 所以當兩者的階碼超過24時相加不會再改變結果. 而0.000061的指數是-15, 因此當res到達pow(2, 10)后0.000061被徹底忽略.