中间浮点精度

# 中间浮点精度 来源：https://randomascii.wordpress.com/2012/03/21/intermediate-floating-point-precision/ 蝙蝠侠，猜猜看：这些计算的精度是多少？ `` void MulDouble(double x, double y, double* pResult) { *pResult = x * y; } void MulFloat(float x, float y, float* pResult) { *pResult = x * y; } `` 如果你回答“double”和“float”，那么你的年轻理想主义得分一分，但正确性得分为零。正确答案，为零理想主义得分并得四十二分正确性得分是“取决于情况”。 阅读下文了解更多细节。 这取决于编译器、编译器版本、编译器设置、32位与64位，以及一些CPU标志的运行时状态。而且，这种“取决于情况”的行为会影响性能，在某些情况下还会影响计算结果。多么令人兴奋！ ## 本系列往期内容…… 如果你刚刚加入，那么阅读本系列前面的文章可能会有所帮助。 - 1:浮点格式技巧 (https://randomascii.wordpress.com/2012/01/11/tricks-with-the-floating-point-format/) – float 格式概述 - 2:愚蠢的浮点技巧 (https://randomascii.wordpress.com/2012/01/23/stupid-float-tricks-2/) – 递增整数表示 - 3:不要用 float 存储 (https://randomascii.wordpress.com/2012/02/13/dont-store-that-in-a-float/) – 关于时间的警示故事 - 3b:它们看起来相等... (https://randomascii.wordpress.com/2012/02/11/they-sure-look-equal/) – 特别奖励帖（不在 altdevblogaday 上） - 4:比较浮点数，2012版 (https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/) - 5:浮点精度——从零到100+位 (https://randomascii.wordpress.com/2012/03/08/float-precisionfrom-zero-to-100-digits-2/) - 5b:C++ 11 std::async 用于快速浮点格式查找 (https://randomascii.wordpress.com/2012/03/11/c-11-stdasync-for-fast-float-format-finding/) – 特别奖励帖（不在 altdevblogaday 上） - 6:中间精度 – (return *this;) 第4帖（比较浮点数）特别相关，因为它的测试代码在不同的编译器上会合法地打印出不同的结果。 ## 默认 x86 代码 最明显的起点是默认的 VC++ 项目。我们来看看 Win32 (x86) 项目 Release 构建的结果，编译设置保持默认，但关闭了链接时代码生成 (LTCG)。以下是生成的代码： > ?MulDouble@@YAXNNPAN@Z push ebp mov ebp, esp fld QWORD PTR \_x$\[ebp\] mov eax, DWORD PTR \_pResult$\[ebp\] fmul QWORD PTR \_y$\[ebp\] fstp QWORD PTR \[eax\] pop ebp ret 0 ?MulFloat@@YAXMMPAM@Z push ebp mov ebp, esp fld DWORD PTR \_x$\[ebp\] mov eax, DWORD PTR \_pResult$\[ebp\] fmul DWORD PTR \_y$\[ebp\] fstp DWORD PTR \[eax\] pop ebp ret 0 有趣的是，MulDouble 和 MulFloat 的代码除了加载和存储指令上的大小说明符外完全相同。这些大小说明符只是指示从内存中加载/存储 float 还是 double。无论哪种方式，fmul 都是以“寄存器精度”执行的。 ## 寄存器精度 默认情况下，32位 VC++ x86 项目会生成 x87 (http://en.wikipedia.org/wiki/X87) 代码来处理浮点运算，如上所示。特殊的 x87 FPU 有八个寄存器，人们通常会告诉你这个 FPU 上的“寄存器精度”意味着 80 位精度。这些人错了。大多数情况下。 事实证明，x87 FPU 有一个精度设置。它可以设置为 24 位 (float)、53 位 (double) 或 64 位 (long double)。VC++ 已经有一段时间不支持 long double 了，因此它在线程启动时将 FPU 初始化为双精度。这意味着，默认情况下，x87 FPU 上的每个浮点运算——加法、乘法、平方根等——都是以双精度完成的，尽管是在 80 位寄存器中。寄存器是 80 位，但每个操作在存储到寄存器之前都会四舍五入到双精度。 但即使这样也不完全正确。尾数被四舍五入到与双精度兼容的 53 位，但指数没有，因此精度与双精度兼容，但范围不同。 在下方的 x87 寄存器图中，蓝色是符号位，粉色是指数，绿色是尾数。始终使用完整的指数，但当舍入设置为 24 位时，只使用浅绿色尾数。当舍入设置为 53 位时，只使用浅绿色和中绿色尾数位，当请求 64 位精度时，使用整个尾数。 80位 x87 寄存器示意图 (https://randomascii.wordpress.com/wp-content/uploads/2012/03/image2.png) 总之，只要你的结果在 DBL_MIN 和 DBL_MAX 之间（应该如此），那么较大的指数就不重要了，我们可以简化一下，只说 x87 上的寄存器精度意味着双精度。 除非不是这样。 虽然 VC++ CRT 将 x87 FPU 的寄存器精度初始化为 _PC_53（53位尾数），但这可以更改。通过调用 _controlfp_s (http://msdn.microsoft.com/en-us/library/e9b52ceh(v=vs.100).aspx) 可以将寄存器值的精度提高到 _PC_64 或降低到 _PC_24。而且，D3D9 有一种习惯，除非你指定 D3DCREATE_FPU_PRESERVE (http://msdn.microsoft.com/en-us/library/windows/desktop/bb172527(v=vs.85).aspx)，否则会将 x87 FPU 设置为 _PC_24 精度。这意味着，在你初始化 D3D9 的线程上，你应该期望许多计算的结果与其他线程不同。D3D9 这样做是为了提高 x87 的浮点除法和平方根的性能。其他操作不会运行得更快或更慢。 因此，上面的汇编代码以以下精度进行中间计算： - 64位精度（80位 long double），如果你以某种方式绕过了 CRT 的初始化，或者如果你使用 _controlfp_s 将精度设置为 _PC_64 - 或者默认情况下为 53位精度（64位 double） - 或者 24位精度（32位 float），如果你初始化了 D3D9 或使用 _controlfp_s 将精度设置为 _PC_24 在更复杂的计算中，编译器可以通过将临时结果溢出到内存来进一步混淆，这可能会以与当前选定的寄存器精度不同的精度进行。 最棒的是，你无法通过查看代码来判断。精度标志是一个每线程的运行时设置，因此同一进程的不同线程调用相同的代码可能会给出不同的结果。这是令人困惑的棒，还是棒到令人困惑？ ## /arch:SSE/SSE2 SSE 的缩略词汤改变了这一切。SSE 家族的指令都指定了使用何种精度，并且没有精度覆盖标志，因此你可以确切知道你将得到什么（只要没有人更改舍入标志，但让我们假装我没提这回事吧？） 如果我们启用 /arch:sse2，其他编译器设置保持不变（release 构建，关闭链接时代码生成，/fp:strict），我们将在 x86 构建中看到这些结果： > ?MulDouble@@YAXNNPAN@Z push ebp mov ebp, esp movsd xmm0, QWORD PTR \_x$\[ebp\] mov eax, DWORD PTR \_pResult$\[ebp\] mulsd xmm0, QWORD PTR \_y$\[ebp\] movsd QWORD PTR \[eax\], xmm0 pop ebp ret 0 ?MulFloat@@YAXMMPAM@Z push ebp mov ebp, esp movss xmm0, DWORD PTR \_x$\[ebp\] movss xmm1, DWORD PTR \_y$\[ebp\] mov eax, DWORD PTR \_pResult$\[ebp\] cvtps2pd xmm0, xmm0 cvtps2pd xmm1, xmm1 mulsd xmm0, xmm1 cvtpd2ps xmm0, xmm0 movss DWORD PTR \[eax\], xmm0 pop ebp ret 0 MulDouble 代码看起来很像，只是三条指令的助记符和寄存器发生了变化。很简单。 MulFloat 代码看起来……更长。更奇怪。更慢。 MulFloat 代码多了四条指令。其中两条将输入从 float 转换为 double，还有一条将结果从 double 转换回 float。第四条是显式加载 'y'，因为当组合 float 和 double 精度时，SSE 无法在一条指令中完成加载和乘法。与 x87 指令集不同，x87 的转换是加载/存储过程的一部分，而 SSE 的转换必须显式进行。这带来了更大的控制权，但增加了成本。如果我们乐观地假设 'x' 和 'y' 的加载和 float 到 double 的转换并行进行，那么这两个函数的浮点运算依赖链差异很大： MulDouble： > movsd -> mulsd -> movsd MulFloat > movss -> cvtps2pd -> mulsd -> cvtpd2ps -> movss 这意味着 MulFloat 的依赖链比 MulDouble 长 66%。movss 和 movsd 指令成本在便宜到免费之间，而转换指令的成本与浮点乘法相当，因此实际延迟增加可能更高。测量这类事情最好情况下也很棘手，而且测量结果外推效果不佳，但在我的粗略测试中，我发现在 VC++ 2010 优化后的 32 位 /arch:SSE2 构建中，MulFloat 的运行时间比 MulDouble 长 35%。在函数调用开销较小的测试中，我看到 float 代码比 double 代码长 78%。呜呼。 ## 扩宽还是不扩宽 宽加载符号 (https://randomascii.wordpress.com/wp-content/uploads/2012/03/image3.png) 那么编译器为什么要这样做？为什么浪费三条指令来将计算扩宽到双精度？ 这样做合法吗？或许是必须的吗？ 关于是否应该使用双精度中间结果进行 float 计算的指导方针多年来一直来回摇摆，除非当它既不向前也不向后的时候。 IEEE 754 浮点数学标准选择不对此做出规定： > 附录 B.1：匿名目标格式由语言表达式求值规则定义。 好的，所以这取决于语言——让我们看看它们怎么说。 C 程序设计语言封面 (http://upload.wikimedia.org/wikipedia/commons/5/54/The_C_Programming_Language_1st_edition_cover.jpg) 在 *The C Programming Language* 中，版权 © 1978（早于 1985 年 IEEE 754 浮点标准），Kernighan 和 Ritchie 写道： > C 语言中的所有浮点运算都以双精度进行；每当 float 出现在表达式中时，它都会通过零填充其小数部分来扩展为 double。 但是 C++ 1998 标准似乎没有包含这段描述。最相关的规定是通常算术转换 (http://msdn.microsoft.com/en-us/library/3t4w2bkb.aspx)（1998 标准中第 5.9 节），它基本上说仅当表达式中包含 double 时才使用双精度： > 5.9：如果上述条件不满足且任一操作数为 double 类型，则另一操作数转换为 double 类型。 但 C++ 1998 标准似乎也允许浮点数学运算以更高精度进行，如果编译器愿意的话： > 5.10 浮点操作数的值和浮点表达式的结果可以以比类型要求的更大的精度和范围来表示；类型不会因此改变。55） C99 标准也似乎允许但不要求 (http://grouper.ieee.org/groups/754/meeting-materials/2001-07-18-c99.pdf) 浮点数学运算以更高精度进行： > 求值类型可能比语义类型更宽 – 宽求值不会拓宽语义类型 而 Intel 的编译器允许你选择中间结果应使用源数的精度、双精度还是扩展精度 (http://cache-www.intel.com/cd/00/00/34/76/347605_347605.pdf)。 > /fp:double: 将中间结果四舍五入到 53 位（double）精度并启用值安全优化。 显然编译器 *可以* 为中间结果使用更高精度，但是 *应该* 吗？ 经典的 Goldberg 文章 (http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) 指出了使用双精度进行 float 计算可能出现的错误： > 这表明，以最高可用精度计算每个表达式并不是一个好规则。 另一方面，微软的 Eric Fleegal 的一篇论证充分文章 (http://msdn.microsoft.com/en-us/library/aa289157%28VS.71%29.aspx) 说： > 通常更可取的做法是尽可能以高精度计算中间结果。 共识是明确的：有时高精度中间结果很棒，有时则很糟糕。IEEE 数学标准参考 C++ 标准，而 C++ 标准参考编译器编写者，编译器编写者有时又参考开发者。 ## 这是一个设计决策和一个性能错误 地板蜡 (https://randomascii.wordpress.com/wp-content/uploads/2012/03/image4.png) 正如 Goldberg 和 Fleegal 所指出的，没有简单的答案。有时更高精度的中间结果会保留重要信息，有时会导致令人困惑的差异。 我想我已经逆向工程了 VC++ 团队如何决定在 /arch:SSE2 下使用双精度中间结果。多年来，Windows 上的 32 位代码一直使用 x87 FPU，并设置为双精度，因此多年来中间值（只要它们留在寄存器中）一直是双精度。显然，当 SSE 和 SSE2 出现时，VC++ 团队希望保持一致性。这意味着显式编码双精度临时变量。 编译器（剧透警告：在 VS 2012 之前）也有一种强烈的倾向，即对任何返回 float 或 double 的函数使用 x87 指令，大概是因为结果必须通过 x87 寄存器返回，而且我确信他们希望避免同一程序内的不一致。 上面我的测试函数中令人沮丧的一点是这些额外指令显然是不需要的。它们对结果不会产生影响。 我如此确信在上述例子中双精度临时变量毫无影响的原因是只有一个操作正在执行。IEEE 标准一直保证基本操作（加、减、乘、除、平方根等）给出完美的舍入结果。在任何单个基本操作上，如果两个 float 的结果存储到 float 中，那么将计算拓宽到双精度是完全无意义的，因为结果已经 *完美* 了。优化器应该识别这一点并移除那四条多余的指令。中间精度在复杂计算中很重要，但在单个乘法中，它不重要。 简而言之，尽管 VC++ 在此情况下的策略是对 float 计算使用双精度，但“as-if”规则意味着编译器可以（并且应该）在速度更快且结果与使用双精度“as-if”相同的情况下使用单精度。 ## 64位改变一切 这是我们的测试代码： `` void MulDouble(double x, double y, double* pResult) { *pResult = x * y; } void MulFloat(float x, float y, float* pResult) { *pResult = x * y; } `` 这是 VC++ 2010 中默认 Release 配置（禁用 LTCG）下为我们的测试函数生成的 x64 机器码： > ?MulDouble@@YAXNNPEAN@Z mulsd xmm0, xmm1 movsdx QWORD PTR \[r8\], xmm0 ret 0 ?MulFloat@@YAXMMPEAM@Z mulss xmm0, xmm1 movss DWORD PTR \[r8\], xmm0 ret 0 差异几乎是喜剧性的。32 位下最短的函数有八条指令，而这些都是三条指令。发生了什么？ 主要区别来自 64 位 ABI (http://blogs.msdn.com/b/freik/archive/2005/03/17/398200.aspx)。在 32 位上，我们需要“push ebp/mov ebp,esp/pop ebp”来设置和拆除堆栈帧以进行快速堆栈遍历。这对于 ETW/xperf 性能分析 (https://randomascii.wordpress.com/category/xperf/) 以及许多其他工具至关重要，并且随着戏剧性的