gcc -mno-sse2 舍入答案

【问题标题】：gcc -mno-sse2 roundinggcc -mno-sse2 舍入
【发布时间】：2016-05-06 06:47:36
【问题描述】：

我正在做一个项目，我将 RGB 转换为 luma，但 -mno-sse2 标志存在一些舍入问题：

这是测试代码：

#include <stdio.h>
#include <stdint.h>

static double rec709_luma_coeff[3] = {0.2126, 0.7152, 0.0722};

int main()
{
    uint16_t n = 242 * rec709_luma_coeff[0] + 242 * rec709_luma_coeff[1] + 242 * rec709_luma_coeff[2];

    printf("%u\n", n);
    return 0;
}

这就是我得到的：

user@gentoo>gcc -mno-sse2 test.c -o test && ./test
241
user@gentoo> gcc test.c -o test && ./test
242

我想 gcc 对 double 乘法使用 sse2 优化，但我不明白为什么优化后的版本是正确的。

另外，您建议我使用什么来获得更一致的结果，ceil() 或 floor()？

【问题讨论】：

与优化无关。无 SSE2 意味着使用旧的 x87 FPU，它比 SSE2 更宽。从某种意义上说，x87 结果的精度更高，但结果可能与使用 SSE2 的结果不同
我建议round() 或nearbyint() 而不是ceil() 或floor()。前者的语义更可能是您想要的。此外，所有这些都可能在其不连续点附近出现一些不稳定，但对于 round() 和 nearbyint()，它们以半整数出现，而对于 ceil() 和 floor()，它们以整数出现。
@user3618511 为什么你首先使用双打进行色彩空间转换？这是史诗般的矫枉过正。
@user3618511 float 也太过分了
@user3618511 这种事情几乎总是用定点数学来完成。例如，luma = (2126 * r + 7152 * g + 722 * b + 5000) / 10000。如果有什么更准确的（你可以让它与浮点数一起工作，但你实际上需要一些带有舍入偏差的黑客）。它也可以用 16 位二进制定点算法进行合理逼近，即使对于标量代码也更有效，并且无限容易与 SIMD 一起使用。

标签： c gcc compilation rounding sse2

【解决方案1】：

TL:DR 使用lrint(x) 或(int)rint(x) 将浮点数转换为整数，并使用舍入到最近而不是截断。不幸的是，并非所有编译器都有效地内联相同的数学函数。见round() for float in C++

gcc -mno-sse2 必须将 x87 用于 double，即使是 64 位代码。 x87 寄存器的内部精度为 80 位，但 SSE2 在 XMM 寄存器中本机使用 IEEE binary64 (aka double) 格式，因此所有临时寄存器在每一步都舍入为 64 位 double。

问题并不像the double rounding problem（80 位 -> 64 位，然后是整数）那么有趣。它也不是来自gcc -O0（默认值：没有额外的优化）在将临时对象存储到内存时进行舍入，因为您在一个 C 语句中完成了整个操作，因此它只对整个表达式使用 x87 寄存器。

只是 80 位精度导致结果略低于 242.0 并被 C 的 float->int 语义截断为 241，而 SSE2 产生略高于 242.0 的结果，截断为 242 . 对于 x87，对于从 1 到 65535 的任何输入，向下舍入到下一个较小的整数是一致的，而不仅仅是 242。（我使用 atoi(argv[1]) 制作了您的程序版本，因此我可以测试其他值，并使用 -O3 )。

记住int foo = 123.99999 是123，因为C 使用“截断”舍入模式（趋向于零）。对于非负数，这与 floor 相同（向 -Infinity 舍入）。 https://en.wikipedia.org/wiki/Floating-point_arithmetic#Rounding_modes.

double 不能准确地表示系数：我用gdb 打印它们并得到：{0.21260000000000001, 0.71519999999999995, 0.0722}。这些十进制表示可能不是以 2 为底的浮点值的精确表示。但它们足够接近，可以看到系数加起来为0.99999999999999996（使用任意精度计算器）。

我们得到一致向下舍入，因为 x87 内部精度高于系数的精度，所以 n * rec709_luma_coeff[0] 等中的和舍入误差，在总结结果时，是 ~2^11 小于差异系数之和与 1.0 之间。（64 位有效位与 53 位）。

真正的问题是 SSE2 版本是如何工作的！据推测，在足够多的情况下，至少在 242 个情况下，临时对象上的舍入到最近偶数恰好向上。它恰好为更多情况产生原始输入，但它为 5、7、10、13、 14, 20...（1..1000 中的前 1000 个数字中有 252 个被 SSE2 版本“删除”，所以它也不总是有效。）

使用-O3 作为您的源代码，它会在编译时以更高的精度进行计算并产生准确的结果。即它与printf("%u\n", n); 编译相同。

顺便说一句，您应该使用static const 作为常量，以便 gcc 可以更好地优化。不过，static 比普通全局要好得多，因为编译器可以看到编译单元中没有任何内容写入值或将它们的地址传递到任何地方，因此它可以将它们视为 const。

【讨论】：