不同编译器的计算结果为何会不一致

我们都知道浮点数存在大数吃小数问题且浮点数的运算不满足结合率。自从 IEEE754 一统江湖之后,理论上指定浮点数数据类型为 double 或者 float 后,同一份 C 代码在不同编译器下的运行结果应该是一样的,但实际并非如此,主要有两方面的原因:

  1. 编译器过于激进的优化,例如 -ffast-math,以及 openmp 等并行方案,他们都忽略了浮点数不满足结合律的特征
  2. 编译器生成的二进制代码中浮点数使用了 X87 或 sse 系列指令集。其中 X87 指令集的中寄存器长度为 80bit = 1 + 15 + 1 + 63,而 sse 符合 IEEE754 规范 64bit = 1 + 13 + 52

经典例子

#include <float.h>
#include <stdio.h>

double maybeTheSame(double x)
{
	return (x + 1.0) - 1.0;
}

float mulSub1(float x, float y, float z)
{
	return x * y - z;
}
float mulSub2(float x, float y, float z)
{
	float xy = x * y;
	return xy - z;
}


int main()
{
#ifdef _MSC_VER
    // https://learn.microsoft.com/zh-cn/cpp/build/reference/arch-x86?view=msvc-170
    // https://learn.microsoft.com/zh-cn/cpp/c-runtime-library/reference/control87-controlfp-control87-2?view=msvc-170
	unsigned int control_word_x87 = 0;
	int result = __control87_2(_PC_64, MCW_PC, &control_word_x87, 0);
#endif

	{
		double a = 1.0 / (1 << 30) / (1 << 23);  // 2^{-53} = DBL_EPSILON / 2 在 IEEE754 下和 1 相加会被大数吃小数
		double r = maybeTheSame(a);
		printf("maybeTheSame(%e) = %e\n", a, r);
	}
	{
		double a = 1.0 / (1 << 30) / (1 << 30) / (1 << 4);  // 2^{-64} 在 X87 指令集下和 1 相加会被大数吃小数
		double r = maybeTheSame(a);
		printf("maybeTheSame(%e) = %e\n", a, r);
	}
	{
		// https://stackoverflow.com/questions/3206101/extended-80-bit-double-floating-point-in-x87-not-sse2-we-dont-miss-it/32699271#32699271
		double a = 1e16;
		double b = 2.9999;
		double c = a + b;
		printf("1e16 + 2.9999 = %.1f\n", c);
	}
	{
		// https://stackoverflow.com/questions/22710272/difference-in-floating-point-arithmetics-between-x86-and-x64
		float a = 50.0f;
		float b = 1.3f;
		float c = 65.0f;
		float r1 = mulSub1(a, b, c);
		printf("mulSub1(%e, %e, %e) = %e\n", a, b, c, r1);
		float r2 = mulSub2(a, b, c);
		printf("mulSub2(%e, %e, %e) = %e\n", a, b, c, r2);
	}

	return 0;
}

上述代码有一些经典的例子。

gcc 测试结果

可以看出

  1. -ffast-math 会忽略结合律,函数 maybeTheSame 直接返回输入的值,因此无论输入 a 多小都返回 a 本身
  2. 在默认设置下,gcc 编译 32 位程序会使用 x87 指令集,编译 64 位程序则使用 sse 指令集。那么我们也可以通过 -mfpmath 倒反天罡,如下图

VS 测试结果

由于 VS 在 64 位时不支持使用 X87 指令集,故仅在 x86 环境下测试

以下是 VS 的结果(采用 sse2 指令集)

在开启如下设置后(/arch:IA32 采用 x87 指令集)

在上述设置的默认情况下结果如下

在开启 X87 精度控制后结果如下

  1. 微软十分保守,即使开启 /fp:fast 也不忽略浮点数不满足结合律的特征,即 maybeTheSame 会照常计算。
  2. 微软提供 __control87_2 可对 x87 指令寄存器的精度进行了控制。
  3. 微软将 x87 指令集的默认长度从 80 改为 64,从而使用双精度 double 时满足 IEEE754,故而结果和 sse2 保持一致;但在这种默认设置下,单精度浮点数 float 的临时结果还是用 64 位保存,故而不满足 IEEE754,最终结果和 sse2 不一致(其实都适用 float 也不那么在乎这点精度误差了)。当然可以通过 __control87_2 设置使 float 的结果也和 sse2 一致。

误差放大

atan2(maybeTheSame(2^{-53}), 0) 在 sse2 指令集中变成 atan2(0, 0) 值为 0,但在 x87 指令集下会变成 atan2(2^{-53}, 0) 值为 π2

posted @   cuzperf  阅读(29)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
点击右上角即可分享
微信分享提示