浮点运算潜在的结果不一致问题

昨天阿楠发现了项目中的一个 bug ，是因为浮点运算的前后不一致导致的。明明是完全相同的 C 代码，参数也严格一致，但是计算出了不相同的结果。我对这个现象非常感兴趣，仔细研究了一下成因。

原始代码比较繁杂。在弄清楚原理后，我简化了出问题的代码，重现了这个问题：

static void
foo(float x) {
    float xx = x * 0.01f;
    printf("%d\n", (int)(x * 0.01f));
    printf("%d\n", (int)xx);
}

int
main() {
    foo(2000.0f);
    return 0;
}

使用 gcc 4.9.2 ，强制使用 x87 浮点运算编译运行，你会发现令人诧异的结果。

gcc a.c -mfpmath=387

19
20

前一次的输出是 19 ，后一次是 20 。

这是为什么呢？让我们来看看 gcc 生成的代码，我截取了相关的段落：

    flds    16(%rbp)
    flds    .LC0(%rip)
    fmulp   %st, %st(1)
    fstps   -4(%rbp)          ; 1. x * 0.01f 结果保存到内存中的 float 变量中
    flds    16(%rbp)
    flds    .LC0(%rip)
    fmulp   %st, %st(1)
    fisttpl -20(%rbp)        ; 2. x * 0.01f 结果直接转换为整型
    movl    -20(%rbp), %eax
    movl    %eax, %edx
    leaq    .LC1(%rip), %rcx
    call    printf
    flds    -4(%rbp)                 ; 3. 读出 1. 保存的乘法结果
    fisttpl -20(%rbp)
    movl    -20(%rbp), %eax
    movl    %eax, %edx
    leaq    .LC1(%rip), %rcx
    call    printf

这里我做了三行注释。

首先，0.01 是无法精确表示成 2 进制的，所以 * 0.01 这个操作一定会存在误差。

两次运算都是 x * 0.01f ，虽然按 C 语言的转换规则，表达式中都是 float 时，按 float 精度运算。但这里 gcc 生成的代码并没有严格设置 FPU 的精度控制，在注释 2 这个地方，乘法结果是直接从浮点寄存器转换为整数的。而在注释 1 这个地方，把乘法结果通过 fstps 以低精度形式保存到内存，再在注释 3 的地方 flds 读回。

所以在注释 2 和注释 3 的地方，浮点寄存器 st 内的值其实是有差别的，这导致了 fisttpl 转换为整数后结果不同。