一個浮點數跨平臺產生的問題

感謝網友唐磊（微博@唐磊_name）投稿，本文原文在唐磊的博客上（原文地址），原文分析還不夠好，並且可能對人有誤導，因此，我對原文作了不少修改，並加了Linux下的內容。浮點數是一個很複雜的事情，但願這篇文章有助於你們瞭解浮點數與其相關的C/C++的編譯選項。（注：我沒有Windows 32位以及C#的環境，因此，對於Windows 32位的程序和C#的程序沒有驗證過）

背景就簡單點兒說，最近一個項目C#編寫，涉及浮點運算，前因後果省去，直接看以下代碼。

1 2 3 4

float p3x = 80838.0f; float p2y = -2499.0f; double v321 = p3x * p2y; Console.WriteLine(v321);

很簡單吧，立刻筆算下結果爲-202014162，沒問題，難道C#沒有產生這樣的結果？不可能吧，開啓Visual Studio，copy代碼試試，果真結果是-202014162。就這樣完了麼？顯然沒有！你把編譯時的選項從AnyCPU改爲x64試試~(服務器環境正是64位滴哦！！)結果竟然邊成了-202014160，對沒錯，就是-202014160。有點不相信，再跑兩遍，仍然是-202014160。呃，想通了，由於浮點運算的偏差，-202014160這個結果是合理的。

爲何合理呢？很正常，由於上面的p3x和p2y是兩個float類型，雖然v321是double，但也是兩個float類型計算完後再轉成double的，float的精度原本也只有7位，因此，對於這個上億的數，天然沒有辦法保證精度。

可是爲何修改CPU的type會有不一樣的效果？嗯，咱們再試試C/C++。

1 2 3 4 5 6 7 8 9 10 11 12 13

#include using namespace std;   int main() {      float p3x = 80838.0f;      float p2y = -2499.0f;      double v321 = p3x * p2y;      std::cout.precision(15);      std::cout << v321 << std::endl;        return 0; }

上面這段C++代碼在不一樣的平臺下的結果以下：

• Windows 32/64位下：-202014160
• Linux 64位下（CentOS 6 gcc 4.4.7）-202014160，
• Linux 32位下（Ubuntu 12.04+ gcc 4.6.3）是：-202014162

合理的結果應該是-202014160，正確的運算結果是-202014162，合理性是浮點精度不夠形成的（文後解釋了合理性）。如果用兩個double相乘可得正確且合理的運算結果（注：把上面C++的程序中的p3x和p2y的類型聲明成double，就能獲得正確的結果，由於double是雙精度的，float是單精度，因此double有足夠的位數存放更多的數位）。可是咱們有點不明白，爲何Linux 32位下，竟然能算出「正確」的數，而不是「合理」的數。

與C++同樣，C#在32位和64位（DEBUG下，這個後面會說）下沒有獲得一致的結果，那咱們來看一下C++/C#的彙編代碼（使用gdb的disassemble /m main 命令，另外下面只顯示 float * float 而後轉成double的那一行代碼的彙編）

Linux平臺下用G++編譯

1 2 3 4 5 6 7

//C++ 32位系統下 Ubuntu 12.04 8       double v321 = p3x * p2y;     0x0804860f <+27>:  flds   0x18(%esp)     0x08048613 <+31>:  fmuls  0x1c(%esp)     0x08048617 <+35>:  fstpl  0x10(%esp)   .......

1 2 3 4 5 6 7

//C++ 64位系統下 CentOS 6 9           double v321 = p3x * p2y;     0x000000000040083c <+24>:    movss  -0x20(%rbp),%xmm0     0x0000000000400841 <+29>:    mulss  -0x1c(%rbp),%xmm0     0x0000000000400846 <+34>:    unpcklps %xmm0,%xmm0     0x0000000000400849 <+37>:    cvtps2pd %xmm0,%xmm0     0x000000000040084c <+40>:    movsd  %xmm0,-0x18(%rbp)

Windows平臺下用Visual Studio編譯

1 2 3 4 5

//C# AnyCPU編譯，Windows VS2012 double v321 = p3x * p2y; 00000049  fld         dword ptr [ebp-40h] 0000004c  fmul        dword ptr [ebp-44h] 0000004f  fstp        qword ptr [ebp-4Ch]

1 2 3 4 5 6

//C# X64位編譯 Windows7 VS2012 double v321 = p3x * p2y;</pre> 009B43B8 movss xmm0,dword ptr [p3x] 009B43BD mulss xmm0,dword ptr [p2y] 009B43C2 cvtss2sd xmm0,xmm0 009B43C6 movsd mmword ptr [v321],xmm0

從上面的彙編代碼能夠看出，不管是Linux和Windows，C++或C# 32位和64對浮點數的彙編指令並不同。 32位生成代碼用的指令是fld/fmul/fstp等，而64位下的使用了movss/mulss/movsd/的指令。看下來，彷佛這個事情和平臺有關係。

咱們繼續調查，咱們發現，其中fld/fmul/fstp等指令是由FPU(float point unit)浮點運算處理器作的，準確的說，是FPU x87指令，FPU在進行浮點運算時，用了80位的寄存器作相關浮點運算，而後再根據是float/double截取成32位或64位，FPU默認上會盡可能減小因爲須要四捨五入帶來的精度問題。可參看浮點運算標準IEEE-754 推薦標準實現者提供浮點可擴展精度格式(Extended precision)，Intel x86處理器有FPU(float point unit)浮點運算處理器支持這種擴展。

非FPU的狀況是用了SSE中128位寄存器(float實際只用了其中的32位，計算時也是以32位計算的)，這就是致使上述問題產生的最終緣由。詳細分析見文末說明。

知道了這一點，咱們能夠man g++ 看一下文檔，咱們能夠找到一個編譯選項叫：-mfpmath，在32位下，這個編譯選項的默認值是：387，也就是x87 FPU指令，在64位下，這個編譯選項的值是sse，也就是使用SSE的指令。因此，就這篇文章中的這個例子而言，若是你在64bits下加上如 -mfpmath=387，你會獲得「正確的」結果，而不是「合理的」結果。

而在VS2012中C++，編譯選項能夠設置(代碼生成中)可選，/fp:[precise | fast | strict]，本例中Release 32位下用precise 或者 strict將獲得合理的結果(-202014160)，fast將產生正確的結果(-202014162), fast debug/release下結果也不同哦(release下才優化了)。64系統下各個結果能夠你們本身去測試下(Debug/Release)，分別看看VS編譯後產生的中間代碼長什麼樣。（陳皓注：個人VS2012在debug編譯下，不管你怎麼設置/fp的參數值，彙編都是同樣的，使用SSE指令，而Release就不同了，可是個人release下看代碼的彙編很是怪異和源代碼對上號，多年不用Windows開發了，對VS的使用僅停留在VC6++/VC2005上）

因此，咱們在從x87 FPU指令向SSE指令作代碼移植的時候，咱們可能會遇到向這樣的浮點數的精度問題，這個精度問題會屢次科學計算中會更糟糕。這個問題並不簡單的只是在32位和64位中的系統出算，這個問題主要仍是看語言編譯器的實現。在更爲高級的語言中，如：C99或Fortran 2003中，引入了「long double」來作可擴展雙精度（Extension Double），這樣就能夠消除更多的精度問題。

下面咱們把程序改爲long double，（注：其中的類型變成long double）

1 2 3 4 5 6 7 8 9 10 11 12 13

#include using namespace std;   int main() {      long double p3x = 80838.0;      long double p2y = -2499.0;      long double v321 = p3x * p2y;      std::cout.precision(15);      std::cout << v321 << std::endl;        return 0; }

用gdb的disassemble /m main你會看到其中的運算的彙編以下（使用了fmlp指令）：

1 2 3 4 5 6

//linux 32位系統 8       long double v321 = p3x * p2y;     0x08048633 <+63>:  fldt   0x10(%esp)     0x08048637 <+67>:  fldt   0x20(%esp)     0x0804863b <+71>:  fmulp  %st,%st(1)     0x0804863d <+73>:  fstpt  0x30(%esp)

1 2 3 4 5 6

//linux 64位系統 8           long double v321 = p3x * p2y;     0x0000000000400818 <+52>:    fldt   -0x30(%rbp)     0x000000000040081b <+55>:    fldt   -0x20(%rbp)     0x000000000040081e <+58>:    fmulp  %st,%st(1)     0x0000000000400820 <+60>:    fstpt  -0x10(%rbp)

咱們能夠看到，32位系統和64位系統使用了一樣的彙編指令（固然，我沒有那麼多物理機，我只是在VMWare Play的虛擬機上測試的，因此上面的示例並不必定適用於全部的地方，另外，C/C++語言和編譯器和平臺有很是大的關係） ，緣由天然是咱們用到了long double這個擴展雙精度的數據類型。（注：若是你用double或float，在Linux上，32位用x87 FPU 指令編譯，而64位用SSE指令編譯）

好了，咱們再回到C#上來，C#的浮點是支持該標準的，其中其官方文檔也提到了浮點運算可能會產生比返回類型更高精度的值（正如上面的返回值精度就超過了float的精度），並說明若是硬件支持可擴展浮點精度的話，那麼全部的浮點運算都將用此精度進行以提升效率，舉個例子x*y/z, x*y的值可能都在double的能力範圍以外了，但真實狀況可能除以z後又能把結果拉回到double範圍內，這樣的話，用了FPU的結果就會獲得一個準確的double值，而非FPU的就是無窮大之類的了。

因此，對於C#來講，你顯然沒法找到一個像C/C++同樣的利用編譯器選項的來解決這個問題的「解決方案」（其實，用編譯器參數是一個僞解決方案）。

並且，要解決這個問題也不是要修改編譯器選項，由於這個問題明顯不是FPU或是SSE的問題，FPU是個過期的技術，SSE纔是合理的技術，因此，若是你不想你的浮點數在計算上有什麼問題，並且你須要精度準確，正確的解決方案不是搞編譯參數，而是——你必定要使用精度更高字節數更多的數據類型，好比：double 或是long double。

另外，你們在寫代碼的時候得保證明際運行環境/測試環境/開發環境的一致性(包括OS架構啊、編譯選項等)啊（尤爲是C/C++ 並且，編譯器上的參數可能會有不少坑，並且有些坑可能會掩蓋你程序中的問題），否則莫名其妙的問題會產生（本文就是開發環境與運行環境不一致致使的問題，糾結了很久才發現是這個緣由）；遇到涉及浮點運算的時候別忘了有多是這個緣由產生的；float/double混用的狀況得特別注意。

Reference：

[1] C# Language Specification Floating point types
[2] Are floating-point numbers consistent in C#? Can they be?
[3] The FPU Instruction Set

附錄

80838.0f * -2499.0f = -202014160.0浮點運算過程的說明

32位浮點數在計算機中的表示方式爲：1位符號位(s)-8位指數位(E)-23位有效數字(M)。
32位Float = (-1)^s * (1+m) * 2^(e-127), 其中e是實際轉換成1.xxxxx*2^e的指數,m是前面的xxxxx(節約1位)

80838.0f = 1 0011 1011 1100 0110.0= 1.00111011110001100*2^16
有效位M = 0011 1011 1100 0110 0000 000
指數位E = 16 + 127 = 143 =  10001111
內部表示 80838.0 =  0 [1000 1111] [0011 1011 1100 0110 0000 000]
= 0100 0111 1001 1101 1110 0011 0000 0000
= 47 9d e3 00 //實際調試時看到的內存值 多是00 e3 9d 47是由於調試環境用了小端表示法法：低位字節排內存低地址端，高位排內存高地址

-2499.0 = -100111000011.0 = -1.001110000110 * 2^11
有效位M = 0011 1000 0110 0000 0000 000
指數位E = 11+127=138= 10001010
符號位s = 1
內部表示-2499.0 = 1 [10001010] [0011 1000 0110 0000 0000 000]
=1100 0101 0001 1100 0011 0000 0000 0000
=c5 1c 30 00

80838.0 * -2499.0 = ?

首先是指數 e = 11+16 = 27
指數位E = e + 127 = 154 = 10011010
有效位相乘結果爲 1.1000 0001 0100 1111 1011 1010 01 //能夠本身動手實際算下
實際中只能有23位，後面的被截斷即1000 0001 0100 1111 1011 1010 01
相乘結果內部表示=1[10011010][1000 0001 0100 1111 1011 101]
= 1100 1101 0100 0000 1010 0111 1101 1101
= cd 40 a7 dd

結果 =  -1.1000 0001 0100 1111 1011 101 *2^27
=  -11000 0001 0100 1111 1011 1010000
=  -202014160
再轉成double後仍是-202014160.

若是是FPU的話，上面的有效位結果不會被截斷，即
FPU結果 = -1.1000 0001 0100 1111 1011 101001 *2^27
= -11000 0001 0100 1111 1011 1010010
= -202014162

全文完，若本文有紕漏之處歡迎指正。