GCC和VC处理浮点数不一致的问题

概述

本文阐述了在工程中遇到的一个精度方法控制不同而导致的问题。

背景介绍

我们看这一段代码：

#include <stdio.h>

int main(int argn, char *argv[])
{
// 上面都是废话，这里正文开始
    double r1, r2;
    int x;

    r1 = 40;
    r2 = 0.85;
    x = (int) (r1 * r2);
    printf("x = %d\n", x);
}

如果用VC编译则输出34，如果用gcc for linux编译则输出33（用gcc for Mac OS X编译也是34）。

完整的分析

为了便于阅读，我这里只给出最相关的资料，如果需要进一步了解需要做试验、检索更详细的资料。

x87 FPU的资料

我们先了解一下FPU的数据寄存器（Data Registers）和控制字（Control Word）。
数据寄存器是80位，有效数字64位。它可以选择几种工作模式：

Precision	PC field	备注
单精度（24位有效数字）	00	32位，即float
保留	01
双精度（53位有效数字）	10	64位，即double
扩展双精度（64位有效数字）	11	80位

PC字段是控制字的8、9两位。
另外，控制字的10、11两位为RC字段，如果为11则表示采用向下取整（Truncate）的方式，而为00则表示四舍五入（Round to near）方式。
注：关于PC、RC的说明，见Book1的8-11页。

计算过程

上述代码的计算过程如下（这是gdb反汇编gcc编译的结果）：

fldl   r2 ; 将r2放入栈
    fmull  r1 ; 和r1相乘，结果在栈顶
    fnstcw cw ; 取FPU的状态字，保存到临时变量cw中
    mov    cw,%ax ; 将cw复制到寄存器ax中
    mov    $0xc,%ah ; 修改RC字段为11，即采用向下取整；修改PC字段为00，即单精度（不过这点并不影响结果）
    mov    %ax,cw2 ; 将修改后的状态字结果放入另一个临时变量cw2中
    fldcw  cw2 ; 更新FPU的状态字
    fistpl x ; 将栈顶的结果取整保存到结果x中

注：我适当改写了gdb反汇编的代码以便于阅读。
注：gdb反汇编标记方法和Intel标准有所不同，第一个操作字是源，第二个则是目的，正好和Intel反过来，需要注意这一点。
因为二进制浮点数是不能准确表达0.85的，会小一点，所以实际40*0.85乘法结果会略小于34，采用向下取整以后就结果为33。

进一步说明

为了便于试验，我将这个计算过程在VC下重新实现一次，代码如下：

#include <stdio.h>

int main(int argn, char *argv[])
{
// 废话结束，正文开始
    double r1, r2;
    unsigned short cw;
    int x;

    r1 = 0.85;
    r2 = 40;

    __asm
    {
        mov cw, 0x037f;

        fldcw cw; // 设置状态字的PC字段为11，即采用扩展双精度，windows原先默认是10，即双精度

        fld r1; // r1入栈
        fld r2; // r2入栈
        fmulp st(1), st // 相乘，结果在栈顶
        fnstcw cw;
        mov ax, cw;

        mov ah, 0x0c;

        mov cw, ax;

        fldcw cw; // 这一段代码同gcc的结果，设置PC为00即单精度，设置RC为11即向下取整

        fistp x; // 取整，保存到x中
    }

    printf("x = %d\n", x);
    return 0;
}

注：在VC下书写汇编代码和Intel的标准一致，目的操作数在前，源在后。
这个结果就是33。
好，我们要讨论的是结果为什么是33？
我们知道，这是因为计算结果略小于34（但是非常接近，在寄存器中只差2^-64次方以内，即最后一位有效数字），而最后采用的是单精度向下取整，自然就少了1，变成33。然而，我们只要把第一句修改为mov cw, 0x027f，这段代码的结果就正确了。为什么？这是因为FPU内部总是采用高精度（80位）进行计算，如果PC为10即双精度，那么FPU在计算结束以后，需要保存为双精度，此时采用了四舍五入（因为RC为00，四舍五入），那么结果将变成34。
所以，因为gcc for linux编译出来的程序默认是80位精度，在最后取整时简单的采用了向下取整，而VC编译出来的程序则是默认64位的精度，在计算过程中就先进行了四舍五入，最后取整时虽然也作了向下取整，但是已经无碍于最终结果了。

再进一步说明

如果gcc将-msse2开关打开，则可以获得满意的结果。让我们看一下相关代码：

...
    fstpl  r3 ; 将结果保存到r3中
    movsd  r3,%xmm0 ; 将r3放到xmm0中
    cvttsd2si %xmm0,%eax ; 取整，双精度浮点->整数
    mov %eax, x ; 将结果放到x中
    ...

第一步即fstpl r3将栈顶的计算结果（80位）保存到内存变量（64位）中，这步操作采用了四舍五入（gcc默认的RC也是00，即四舍五入）将80位精度转为64位精度，往后自然全部正确。
按照这个道理，我们可以用这样的方法改写C代码获得正确结果：

...
    r3 = r1 * r2;
    x = (int) r3;
    ...

这样结果也正确，因为r3 = r1 * r2这一步会进行四舍五入。

结论

gcc for linux的浮点数直接取整实现是一个bug。

最后

浮点数计算对完美的移植是一个重大的考验，因为对浮点处理方法的不一致性、或处理顺序的不一致性会导致精度上的微小差异。这个微小的差异随着计算的进行，会逐步的放大。

相关文档

Book1：Intel 64 and IA-32 Architectures Software Developer's Manual VOLUME 1: Basic Architecture