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

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

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

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

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

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

这里我做了三行注释。

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

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

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


2018 年 6 月 14 日补充:

在 DDJ 1997 年的一篇访谈文章中,谈及了类似的问题。

引用如下:

Let me give an example, which arose yesterday. A student was doing a computation involving a simulation of a plasma. It turns out that as you go through this computation there will be certain barriers. This may be a reflecting barrier. That means if a particle moves through this barrier, it should not really go through, it should be reflected. Others may be absorbing, others may be periodic. He has a piece of code that is roughly

As far as we can tell, when he turns on optimization, the value of x is computed in a register with extra width. This happens on the Intel machine because the registers have extra width. It also happens on some others. The value of x that he computes is not, in fact, a float. It's in a register that's wider. The machine stores x somewhere, and in the course of doing that, converts it to a float. But the value used in the comparison is the value in the register. That x is wider. The condition that the register x be greater than or equal to j doesn't guarantee that x, when stored in y, will be less than j. Sometimes y=j, and that should never be. My student counted on what compiler writers call "referential transparency" and was disappointed because the compiler writer saved some time in optimization by not reloading the register from the stored value.

云风 提交于 July 25, 2017 10:31 AM | 固定链接

COMMENTS

大神,同样的问题。在C#里,如果Debug编译,选择AnyCPU,结果分别是20、19.如果Release编译,结果就是19、19

Posted by: Yao | (16) October 30, 2017 05:49 PM

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa

Posted by: asd123 | (15) September 18, 2017 07:12 PM

这种浮点型数据运算结果不一致有什么办法规避或者解决呢?

Posted by: ahuang | (14) August 9, 2017 03:35 PM

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323

Posted by: Nothing | (13) July 27, 2017 11:49 PM

确实是与优化选项和编译器版本有很大的关系。弄出代码完全一致,但在不同编译器或者不同编译选项上输出不同结果的浮点运算是很容易的。

Posted by: ghosting | (12) July 27, 2017 08:48 AM

@ghosting

这个问题和 c 代码是否完全一致关系不大。构造出一个c代码完全一致,但结果有差别的例子也是可以的,就是比较繁琐,我懒得做。而且和优化选项,gcc 版本关系更大,很难重现。

Posted by: cloud | (11) July 27, 2017 12:35 AM

两次运算不是完全相同的C代码,浮点数并不是实数,这个结果也没有什么可以诧异的。
C/C++对于浮点数的定义很宽松的,标准上甚至没有定义过float必须是32位。
x87这种硬件实际上是支持Extended precision格式,但不完全支持single/double格式(只支持存储相关的操作)。
VC对这种情况加了个strictfp的扩展,使得每次操作的中间结果都会转换下,gcc没有等价的玩意。
本来写出没有问题的浮点代码就有很多点需要注意,在x87上因为硬件的关系会更加棘手一点,最简单的办法还是用strictfp模式或者直接制定fpu为sse。

Posted by: ghosting | (10) July 26, 2017 08:16 PM

@dwing

这只是一个简化的例子。和在 C 语言中是否明确将中间结果暂存在局部变量里无关。

两个完全相同的长表达式,如果有中间结果相同,编译器一样会临时放在内存中的。

实际的代码中大约是

x * 0.01f / a + b 这样的表达式被 inline 后重复了两次,而 由于连续做相同运算,以及恰巧的上下文关系, x * 0.01f 的计算结果被暂存起来了。代码生成和上下文关系很大,不方便举例罢了。

另外,四舍五入并不能解决这个问题,四舍五入只是换了一个截断点而已,把 [1,2) 取整成 1 和把 [0.5,1.5) 取整成 1 只是在数轴上偏移了 0.5 。一样可以制造出不一致。

Posted by: Cloud | (9) July 26, 2017 03:21 PM

多谢dwing,原来C89/90就已经不提升精度了。

Posted by: stirp | (8) July 26, 2017 10:54 AM

常写浮点运算程序的人对取整操作是很敏感的,一般是能不取整就不取整,实在要取整肯定要想想把连续的值离散化的意义,是否应该四舍五入.
另外追求完全一致的计算本来就不应该用浮点数.

Posted by: dwing | (7) July 26, 2017 10:01 AM

这两次运算不能叫完全相同的C代码吧.
虽然从语言上定义此表达式的精度是float,但实际上计算结果的临时数值是尽可能保留精度的.
x87寄存器是80位,从计算结果转成int实际上是从80位精度转换,而中间用float变量保存再去取整就会有细微差异,换成double也可能不一致.

Posted by: dwing | (6) July 26, 2017 09:55 AM

C 语言规范:表达式里最高精度是 float ,那么表达式的值精度就是 float 。没有相乘提升到 double 的说法。

第一个写成 (int)((float)(x * 0.01f)) 也是一样的,并不是对 double 取整的问题。

Posted by: Cloud | (5) July 25, 2017 09:07 PM

刚刚被小伙伴提示默认情况下float相乘是自动提升精度到double的,因此第一次输出就应该是对double取整的。

x87 浮点运算的效果是float相乘结果还是float么?本地不支持387,没法测试。

Posted by: stirp | (4) July 25, 2017 06:57 PM

gcc 4.8.4没有重现

Posted by: alphachen | (3) July 25, 2017 03:12 PM

这种问题不注意真是很意外呢。还和编译器版本相关,gcc 7.1.1 没有复现。

Posted by: 依云 | (2) July 25, 2017 01:47 PM

所以这种计算还是需要四舍五入,强制转换成整形,就直接截断了,估计就是19.999999跟20.000001的差别

参考链接


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

Parsing C++ in Python with Clang

macOS Mojave 10.14.2 llvm 7.0.1

参考链接


GCC支持在代码中对头文件是否存在的判断(__has_include)

我们在实际编写代码的时候,经常需要判断当前编译环境是否存在我们需要的头文件,如果不存在,则使用其他头文件代替。

以前这个操作都是通过外部的configure文件生成Makefile的时候指定。

最近的GCC已经增加了__has_include这个内置函数判断头文件是否存在。

这个功能最早是Clang实现的,现在GCC终于补上了这个功能。

例子如下:

Code Composer Studio 8.2.0.00007 ,GCC 5支持这个内置函数。

参考链接


x86架构实现ARM架构下的SSAT指令

最近在使用Matlab仿真ARM下的Q15,Q31相关的操作,涉及到翻译CMSIS库中的部分函数翻译到x86下运行的情况。

一般会遇到两个比较特殊的宏,一个是__CLZ宏,另一个是 __SSAT宏,前者直接使用__buildin_clz替换就可以非常正常的工作,后者就比较复杂。


可以使用如下方式来翻译这个宏:

int在32/64位系统中的字节数

 

64-bit data models
Data model short (integer) int long (integer) long long pointers,
size_t
Sample operating systems
SILP64 64 64 64 64 64 Classic UNICOS[41] (versus UNICOS/mp, etc.)
ILP64 16 64 64 64 64 HAL Computer Systems port of Solaris to the SPARC64
LLP64,
IL32P64
16 32 32 64 64 Microsoft Windows (x86-64 and IA-64) using Visual C++; and MinGW
LP64,
I32LP64
16 32 64 64 64 Most Unix and Unix-like systems, e.g., Solaris, Linux, BSD, macOS. Windows when using Cygwin; z/OS

参考链接


32位和64位系统区别及int字节数
64-Bit Programming Models: Why LP64?
64-bit computing
What decides the sizeof an integer?

CMake如何编译CUDA(.cu)源文件

现在的项目,如果需要用到计算加速,NvidiaCUDA往往是首选。那么如何在CMake中编译写好的CUDA源代码,可以参考如下。

首先使用FIND_PACKAGE找到已经安装的CUDA,此时需要配置的环境变量等,应该已经自动配置完成了

接下来,使用CUDA_ADD_LIBRARY取代原来的ADD_LIBRARY,如下:

如果是可执行程序,请使用CUDA_ADD_EXECUTABLE取代ADD_EXECUTABLE

参考链接


CMake: how to add cuda to existing project

在Makefile中检查gcc版本

问题描述

我想使用一些在旧的GCC版本中不可用的GCC警告开关(例如-Wtype-limits)。

有没有一个简单的方法来检查GCC版本?

最佳解决方案

我不会说容易,但是可以使用GNU makeSHELL函数来执行gcc --version命令,然后使用ifeq条件表达式来检查版本号,并适当地设置CFLAGS变量。

以下是一个简单的MakeFile示例:

注意:不存在ifgt语法。但是可以使用expr命令进行比较。

例子如下:

次佳解决方案

要将完整的3部分GCC版本(不仅第一位数字)转换成数字格式,适合比较(例如40701)使用

其中解决了任何版本部分中双数字版本号的可能性,以及gcc -dumpversion输出内容缺少3个部分的可能性(在一些较早的GCC版本中存在这种情况)。

所以要测试MakeFile中的版本,使用(注意$$里面的最后一个sed命令)

第三种解决方案

我刚刚遇到这个问题,我需要测试GCC的前两位数字,想要一个比上面sed过滤字符串更复杂的操作。我使用bc进行比较,因为它支持浮点数(expr将非整数(non-integers)都视为字符串):

如果在GCC 4.9之后发布了GCC 4.10,那么用sed处理是必要的:

参考链接


在Makefile中检查gcc版本?

Ubuntu 16.04系统上Clang与GCC之间切换

在编译C++代码的时候,我们有时需要比较一下不同编译器之间优化性能的差异,因此需要在ClangGCC之间进行切换,用来比较最后的实际效果。

Ubuntu 16.04系统上使用如下命令进行切换

参考链接


Switching between GCC and Clang/LLVM using CMake

Ubuntu 16.04上Linux C++程序性能分析工具perf使用入门

目前,perfLinux系统上最全面最方便的一个性能检测工具。由Linux内核携带并且同步更新。

Ubuntu 16.04系统上需要执行如下命令安装:

使用方法如下(gcc编译时最好使用-g参数,生成符号,方便调试):

参考链接