一个浮点数跨平台产生的问题
感谢网友唐磊（微博@唐磊_name）投稿，本文原文在唐磊的博客上（原文地址），原文分析还不够好，而且可能对人有误导，所以，我对原文做了很多修改，并加了Linux下的内容。浮点数是一个很复杂的事情，希望这篇文章有助于大家了解浮点数与其相关的C/C++的编译选项。（注：我没有Windows 32位以及C#的环境，所以，对于Windows 32位的程序和C#的程序没有验证过）
背景就简单点儿说，最近一个项目C#编写，涉及浮点运算，来龙去脉省去，直接看如下代码。
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++。
#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++编译
//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) .......
//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编译
//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]
//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）
#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指令）：
//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)
//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
全文完，若本文有纰漏之处欢迎指正。
（转载本站文章请注明作者和出处 宝酷 – sou-ip ，请勿用于任何商业用途）
(21 人打了分，平均分： 3.57 )
Loading...
《一个浮点数跨平台产生的问题》的相关评论
当初在学校学习基础理论的时候就知道这个754标准，但是当初对它的理解还只停留在这个数据类型表示的是个近似值，判断相等要使用相减小于一个非常小的数。后来实际工作中很少用到浮点型，基本都是在用decimal这样的十进制类型。原来这里面坑还是很深啊，这样的问题也得对汇编，操作系统，编译原理这些底层的基础知识吃得很透才能找到解决问题的方向。从这个角度来说，大学那些计算机科学基础课程还是非常重要的，只是当初不了解它们的意义，学习的时候也算理解，但不够深入。有时间还是要把这些东西从头再过一遍。所谓温故而知新。
BTW：这些厂商号称的一次编写，多处执行还是一个美好的远望，虽然解决的绝大多数问题，跟实际环境还是有一定的出入，不能迷信他们啊。
在FeedlyReader中看到，特地过来拍砖。本文技术细节很清晰，逻辑也很严密，但结论有误导听众的嫌疑。首先，在Intel的指令集参考手册中明确表示，X87已经是过时的指令，不支持超标量和SIMD，推荐使用SSE代替；其次，该例子中，程序员的错误在于没有理解C语言的隐式类型转换规则，把两个单精度相乘赋给双精度，其过程是两个单精度按单精度乘法计算，得到单精度结果，然后隐式转化为双精度赋给目标地址。这个计算本身是损失精度的，只有7到位有效数字。标准解决方法是使用双精度进行计算，只需要把一个操作数强制转换为双精度即可。最后，文章Linux下64位反汇编用错了代码，已经全是双精度了，汇编也是molsd，SSE双精度乘法。事实上，只要硬件平台符合IEEE754浮点标准，软件按所需精度选择类型，其他的留给编译器，浮点精度就实现跨平台啦，GPU和一些FPGA都支持IEEE754了，X86系列CPU从引入SSE2后也完整支持，至于那个X87，历史遗留产物，由它去吧。
@ProgramFan 是的，这个本来就是程序代码的问题，而不应该用编译器选项来解决。本文是要想说的是，编译器参数中的坑。文中特意增强（加粗飘红）了一下结论，希望本文不会产生更多的误导。
一直规避使用float，由于涉及浮点运算的也少，暂时没出过问题；文章不错，再次学习了大学的知识，赞
赞，终于改过来了。一直订阅你的博客，学到很多东西。CentOS下64位反汇编代码仍然有问题，阅读器里的版本是这样的：
//C++ 64位系统下 CentOS 6
……
6 double p3x = 80838.0f;
0x00000000004007ec : movabs $0x40f3bc6000000000,%rax
0x00000000004007f6 : mov %rax,-0x18(%rbp)
7 double p2y = -2499.0f;
0x00000000004007fa : movabs $0xc0a3860000000000,%rax
0x0000000000400804 : mov %rax,-0x10(%rbp)
8 double v321 = p3x * p2y;
0x0000000000400808 : movsd -0x18(%rbp),%xmm0
0x000000000040080d : mulsd -0x10(%rbp),%xmm0
0x0000000000400812 : movsd %xmm0,-0x8(%rbp)
……
其中 6，7行声明的是 double，第八行的汇编也是双精度操作，这应该是正确的程序，而不是文中出问题的版本，建议修正。
@ProgramFan 谢谢，代码搞乱了。
竟然看到了同学（隔壁班）的名字。哈哈
奉劝楼主少玩windows，误终身的平台。
数据类型选择错了，都去调整compiler 参数起来了。这是南辕北辙的节奏。
不知道耗子叔有没有用过OpenCL，在分别使用CPU和GPU进行浮点数运算时精度损失特别明显。@宝酷
interesting
这里说的改编译参数是希望得到合理的结果(某种意义上讲有误差的结果才是正确的结果，默认情况下编译器自作多情给我了个精确的答案)，而不是说想通过编译参数来提高精度。
@Dimtry 这还挑起操作系统滴。。。你用linux写桌面游戏试试。。。就算你会写 也不见得有人用啊。
另改编译参数是希望得到合理的结果而不是想用它提高精度。
@Hallelujah float/double之类各有各的适用地方啊。自己应该在速度/精度上有一个权衡。
linux 还是很有游戏前途的 steam平台专门出了个游戏的linux发行版 移植了各种游戏过去 还准备打造成 xbox这类的动向 未来还是很有看头@tl3shi
作为高级语言用户，不应该对浮点数的精度有过多假设。如果某个具体实现给我们提供了超出通常水平的浮点精度，我们可以接受它，而不应该依赖它。
得出 float 小数“精度有限”的结论就足够了，再往后的分析就有点走火入魔了（入了偏门）。
博主分析的很好。
博主的结论中如果再加上@ProgramFan 的评论，
这篇博文就完美了。
从代码安全的角度看，这里的问题和
uint8_t test_a = 127;
uint8_t test_b = 130;
uint16_t test_result = test_a + test_b;（在这段代码在Cygwin运行正常）
存在数据截断一样。
根据MISRA-C 规定，对于隐式类型转换要慎用，尽量在同类型间操作。
对于C语言未定义的，不予使用。就像char，不知道是有符号还是无符号的。所以还是用uint8_t或者int8_t。
博主的文章对我很有用。谢谢。
不过我遇到的问题，精度解决不了的。
float fVal = 0.69f;
int a = (int)(fVal * 1000);
在VS2010下，得到的a是689。
double dVal = 0.69;
int a = (int)(dVal*10000000);
得到a是6899999。
我突然想起来这篇文章http://www.agner.org/optimize/
有一些很细节的32位和64位的区别.
@ProgramFan
你说的很对，完全是程序的问题