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

weixin_34162629發表於2018-03-03

感謝網友唐磊（微博@唐磊_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

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

相關文章