一个浮点数跨平台产生的问题

感谢网友唐磊投稿,本文原文在唐磊的博客上(原文地址),我对原文做了很多修改,并加了Linux下的内容。浮点数是一个很复杂的事情,希望这篇文章有助于大家了解浮点数与其相关的C/C++的编译选项。(注:我没有Windows 32位以及C#的环境,所以,对于32位的程序和C#的程序没有验证过)

背景就简单点儿说,最近一个项目C#编写,涉及浮点运算,来龙去脉省去,直接看如下代码。

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

很简单吧,马上笔算下结果为-202014162,没问题,难道C#没有产生这样的结果?不可能吧,开启VisualStudio,copy代码试试,果然结果是-202014162。就这样完了么?显然没有!你把编译时的选项从AnyCPU改成x64试试~(服务器环境正是64位滴哦!!)结果居然边成了-202014160,对没错,就是-202014160。有点不相信,再跑两遍,仍然是-202014160。呃,想通了,因为浮点运算的误差,-202014160这个结果是合理的。嗯,再试试C++。

#include <iostream>
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和p3y的类型声明成double,就能得到正确的结果)。

用C++同样的代码,在32位和64位(DEBUG下,这个后面会说)下没有得到一致的结果,那我们来看一下他们的汇编代码(使用gdb的disassemble /m main 命令)

Linux平台下用G++编译

//C++ 32位系统下 Ubuntu 12.04
......
6	    float p3x = 80838.0f;
   0x080485fd <+9>:	mov    $0x479de300,%eax
   0x08048602 <+14>:	mov    %eax,0x18(%esp)

7	    float p2y = -2499.0f;
   0x08048606 <+18>:	mov    $0xc51c3000,%eax
   0x0804860b <+23>:	mov    %eax,0x1c(%esp)

8	    double v321 = p3x * p2y;
   0x0804860f <+27>:	flds   0x18(%esp)
   0x08048613 <+31>:	fmuls  0x1c(%esp)
   0x08048617 <+35>:	fstpl  0x10(%esp)

.......
//C++ 64位系统下 CentOS 6
......
6           double p3x = 80838.0f;
   0x00000000004007ec <+8>:     movabs $0x40f3bc6000000000,%rax
   0x00000000004007f6 <+18>:    mov    %rax,-0x18(%rbp)

7           double p2y = -2499.0f;
   0x00000000004007fa <+22>:    movabs $0xc0a3860000000000,%rax
   0x0000000000400804 <+32>:    mov    %rax,-0x10(%rbp)

8           double v321 = p3x * p2y;
   0x0000000000400808 <+36>:    movsd  -0x18(%rbp),%xmm0
   0x000000000040080d <+41>:    mulsd  -0x10(%rbp),%xmm0
   0x0000000000400812 <+46>:    movsd  %xmm0,-0x8(%rbp)
......

Windows平台下用Visual Studio编译

//C++ 32位系统,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++ 64位系统 Windows7 VS2012
    long double v321 = p3x * p2y;
008443B8  movss       xmm0,dword ptr [p3x]
008443BD  mulss       xmm0,dword ptr [p2y]
008443C2  cvtss2sd    xmm0,xmm0
008443C6  movsd       mmword ptr [v321],xmm0 

从上面的汇编代码可以看出,无论是Linux和Windows,C++ 32位和64对浮点数的指令并不一样,在相同32位或64位上的汇编基本上是一样的。 32位生成代码用的指令是fld/fmul/fstp等,而64位下的使用了movsd/mulsd/movsd以及Windows的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位计算的),这就是导致上述问题产生的最终原因。详细分析见文末说明。虽然SSE的精度不如FPU,但是SSE可以做并行处理,其性能上要比x87 FPU要好。

知道了这一点,我们可以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 <iostream>
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++语言和编译器和平台有非常大的关系) 

好了,我们再回到C#上来,C#的浮点是支持该标准的,其中其官方文档也提到了浮点运算可能会产生比返回类型更高精度的值(正如上面的返回值精度就超过了float的精度),并说明如果硬件支持可扩展浮点精度的话,那么所有的浮点运算都将用此精度进行以提高效率,举个例子x*y/z, x*y的值可能都在double的能力范围之外了,但真实情况可能除以z后又能把结果拉回到double范围内,这样的话,用了FPU的结果就会得到一个准确的double值,而非FPU的就是无穷大之类的了。

但对于C#来说,目前还没找到解决方案

所以大家在写代码的时候得保证实际运行环境/测试环境/开发环境的一致性(包括OS架构啊、编译选项等)啊,不然莫名其妙的问题会产生(本文就是开发环境与运行环境不一致导致的问题,纠结了好久才发现是这个原因);遇到涉及浮点运算的时候别忘了有可能是这个原因产生的;另外,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

全文完,若本文有纰漏之处欢迎指正。

(转载本站文章请注明作者和出处 酷 壳 – CoolShell.cn ,请勿用于任何商业用途)

——=== 访问 酷壳404页面 寻找遗失儿童。 ===——

相关文章

本条目发布于。属于开发分类。作者是
分享到:

发表评论

电子邮件地址不会被公开。 必填项已用*标注