Binary64 浮点加法舍入模式错误和行为差异 32/64 位
Binary64 floating point addition rounding mode error and behaviors difference 32/64 bits
当我尝试在 Intel 核心 I7 / I5 上添加以下两个浮点数时，我注意到一个舍入错误：
2.500244140625E+00 + 4503599627370496.00
<=> 0x1.4008p+1 + 0x1.0p+52
通过 faddl 汇编指令（当我使用 32 位编译器编译时）使用两个 double 精度常量进行添加。
我得到的结果是：
4.503599627370498E+15 = 0x1.0000000000002p+52
而不是：
4.503599627370499E+15 = 0x1.0000000000003p+52
（如我所料，并得到 http://weitz.de/ieee/ 的确认。）
示范：
0x1.0p+52 = 0x10000000000000.00p+0
0x1.4008p+1 = 0x2.801p+0
0x10000000000000.00p+0 + 0x2.801p+0 = 0x10000000000002.801p+0 (exactly)
0x10000000000002.801p+0 = 0x1.0000000000002801p+52 (exactly)
0x10000000000002.801p+0 = 0x1.0000000000003p+52 (after rounding)
我在调试模式下仔细检查并验证我的 FPU 处于 "round to the nearest mode"。
更奇怪的是当我用64位编译器编译我的代码，然后使用addsd指令时，有no舍入错误.
有没有人可以给我参考或解释关于在相同 FPU 但使用不同指令集的 'double' 加法的精度差异？
首先，您正在查看以 10 为基数的数字。您想讨论浮点数和舍入，这需要以 2 为基数进行讨论。
第二个单数和双数具有不同长度的尾数，因此显然对于相同的数字，您四舍五入的位置在小数点 1.2345678 中有所不同，我们可以将其四舍五入为 1.23 或将其四舍五入为 1.2346，具体取决于我们允许一轮向上舍入的位数向下，遵循向上舍入规则。
由于您在这里的某个时刻是以 10 为基数的，因此您还混合了可能的编译时转换、运行 时间操作和 运行 时间转换
我接受
float x=1.234567;
x=x*2.34;
printf("%f\n",x);
有编译时转换，首先是最基本的 ascii 到 double，然后是 double 到 float，以完全符合语言的要求（没有将 F 放在常量的末尾）。然后 运行 时间相乘，然后 运行time 转换为 ascii，运行time C 库可能与编译时间不同，它们是否遵循相同的舍入设置等。很容易找到您只需声明 x=1.234 的数字...然后下一行代码是 printf 而 printf 不是您输入的内容，除了 运行 之外没有浮点数学时间浮点数。
因此，在您提出这个问题之前，我们需要查看您的数字的二进制版本，您的问题的答案应该几乎会自动从中脱颖而出，而无需进一步的帮助，但如果您仍然需要帮助，那么 post那我们可以看看它。进行基于小数的讨论会增加编译器和库问题，并且在出现问题时更难隔离问题。
FPU 寄存器为 80 位宽，每当单精度或双精度数加载 fld 及其变体时，默认情况下它会转换为 double extended precision1.
因此 fadd 通常适用于 80 位数字。
SSE 寄存器与格式无关，SSE 扩展不支持双扩展精度。
例如，addpd 适用于双精度数字。
默认的舍入模式是舍入到最近的（偶数），这意味着通常的舍入到最近的但是为了以防万一平局（例如 4.5 => 4）。
为了实现 IEEE 754 对无限精度数字执行算术的要求，硬件需要两个保护位和一个粘性位2
双
我会写一个双精度数作为
<sign> <unbiased exponent in decimal> <implicit integer part> <52-bit mantissa> | <guard bits> <sticky bit>
两个数
2.500244140625
4503599627370496
是
+ 1 1 0100000000 0010000000 0000000000 0000000000 0000000000 00
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 00
第一个移位
+ 52 0 0000000000 0000000000 0000000000 0000000000 0000000000 10 |10 1
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 00 |00 0
求和完成
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 10 |10 1
四舍五入到最接近的（偶数）得到
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 11
因为 0 |10 1 比 0 |00 0 更接近 1 |00 0。
双扩展
两个数是
+ 1 1 0100000000 0010000000 0000000000 0000000000 0000000000 0000000000 000
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 000
第一个移位
+ 52 0 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000 | 10 0
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 000 | 00 0
求和完成
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000 | 10 0
四舍五入到最接近的（偶数）：
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000
因为 0 | 10 0 是最接近的偶数。
当此数字随后从扩展双精度转换为双精度时（由于 fstp QWORD []），使用扩展尾数的第 52、53 和 54 位重复舍入作为守卫和粘性位
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 10|100
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 10
因为 0|100 再次被打破到最接近的偶数。
1 请参阅英特尔手册第 1 卷第 8.5.1.2 章。
2 保护位是在移动其中一个数以使指数匹配后保留的额外精度位。粘性位是比最不重要的位更重要的位的或。有关格式方法，请参阅 this page and Goldberg 的 "on Rounding" 部分。
感谢我的问题收到的所有评论，我明白了发生了什么并能够解决问题。
这里我尽量总结一下。
首先，确认四舍五入不正确。如 @MarkDickinson, it can be due to a "double rounding", but I do not know if it can be confirmed. Indeed, it can also be due to others phenomenon such as the ones described in the publication given by Pascal Cuoq.
所述
好像ia32 FPU在某些数字四舍五入的问题上并没有完全符合IEEE754标准。
默认情况下，GCC（32 位版本）生成使用 FPU 计算 Binary64 数字加法的代码。
但是，在我的电脑（intel core i7）上，SSE 单元也可以进行这些计算。 GCC（64 位版本）默认使用此单元。
在 GCC32 命令行上使用以下两个选项解决了我的问题。
-msse2 -mfpmath=sse.
(谢谢EOF)
当我尝试在 Intel 核心 I7 / I5 上添加以下两个浮点数时，我注意到一个舍入错误：
2.500244140625E+00 + 4503599627370496.00 <=> 0x1.4008p+1 + 0x1.0p+52
通过 faddl 汇编指令（当我使用 32 位编译器编译时）使用两个 double 精度常量进行添加。
我得到的结果是：
4.503599627370498E+15 = 0x1.0000000000002p+52
而不是：
4.503599627370499E+15 = 0x1.0000000000003p+52
（如我所料，并得到 http://weitz.de/ieee/ 的确认。）
示范：
0x1.0p+52 = 0x10000000000000.00p+0
0x1.4008p+1 = 0x2.801p+0
0x10000000000000.00p+0 + 0x2.801p+0 = 0x10000000000002.801p+0 (exactly)
0x10000000000002.801p+0 = 0x1.0000000000002801p+52 (exactly)
0x10000000000002.801p+0 = 0x1.0000000000003p+52 (after rounding)
我在调试模式下仔细检查并验证我的 FPU 处于 "round to the nearest mode"。
更奇怪的是当我用64位编译器编译我的代码，然后使用addsd指令时，有no舍入错误.
有没有人可以给我参考或解释关于在相同 FPU 但使用不同指令集的 'double' 加法的精度差异？
首先，您正在查看以 10 为基数的数字。您想讨论浮点数和舍入，这需要以 2 为基数进行讨论。
第二个单数和双数具有不同长度的尾数，因此显然对于相同的数字，您四舍五入的位置在小数点 1.2345678 中有所不同，我们可以将其四舍五入为 1.23 或将其四舍五入为 1.2346，具体取决于我们允许一轮向上舍入的位数向下，遵循向上舍入规则。
由于您在这里的某个时刻是以 10 为基数的，因此您还混合了可能的编译时转换、运行 时间操作和 运行 时间转换
我接受
float x=1.234567;
x=x*2.34;
printf("%f\n",x);
有编译时转换，首先是最基本的 ascii 到 double，然后是 double 到 float，以完全符合语言的要求（没有将 F 放在常量的末尾）。然后 运行 时间相乘，然后 运行time 转换为 ascii，运行time C 库可能与编译时间不同，它们是否遵循相同的舍入设置等。很容易找到您只需声明 x=1.234 的数字...然后下一行代码是 printf 而 printf 不是您输入的内容，除了 运行 之外没有浮点数学时间浮点数。
因此，在您提出这个问题之前，我们需要查看您的数字的二进制版本，您的问题的答案应该几乎会自动从中脱颖而出，而无需进一步的帮助，但如果您仍然需要帮助，那么 post那我们可以看看它。进行基于小数的讨论会增加编译器和库问题，并且在出现问题时更难隔离问题。
FPU 寄存器为 80 位宽，每当单精度或双精度数加载 fld 及其变体时，默认情况下它会转换为 double extended precision1.
因此 fadd 通常适用于 80 位数字。
SSE 寄存器与格式无关，SSE 扩展不支持双扩展精度。
例如，addpd 适用于双精度数字。
默认的舍入模式是舍入到最近的（偶数），这意味着通常的舍入到最近的但是为了以防万一平局（例如 4.5 => 4）。
为了实现 IEEE 754 对无限精度数字执行算术的要求，硬件需要两个保护位和一个粘性位2
双
我会写一个双精度数作为
<sign> <unbiased exponent in decimal> <implicit integer part> <52-bit mantissa> | <guard bits> <sticky bit>
两个数
2.500244140625
4503599627370496
是
+ 1 1 0100000000 0010000000 0000000000 0000000000 0000000000 00
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 00
第一个移位
+ 52 0 0000000000 0000000000 0000000000 0000000000 0000000000 10 |10 1
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 00 |00 0
求和完成
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 10 |10 1
四舍五入到最接近的（偶数）得到
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 11
因为 0 |10 1 比 0 |00 0 更接近 1 |00 0。
双扩展
两个数是
+ 1 1 0100000000 0010000000 0000000000 0000000000 0000000000 0000000000 000
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 000
第一个移位
+ 52 0 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000 | 10 0
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 000 | 00 0
求和完成
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000 | 10 0
四舍五入到最接近的（偶数）：
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000
因为 0 | 10 0 是最接近的偶数。
当此数字随后从扩展双精度转换为双精度时（由于 fstp QWORD []），使用扩展尾数的第 52、53 和 54 位重复舍入作为守卫和粘性位
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 1010000000 000
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 10|100
+ 52 1 0000000000 0000000000 0000000000 0000000000 0000000000 10
因为 0|100 再次被打破到最接近的偶数。
1 请参阅英特尔手册第 1 卷第 8.5.1.2 章。
2 保护位是在移动其中一个数以使指数匹配后保留的额外精度位。粘性位是比最不重要的位更重要的位的或。有关格式方法，请参阅 this page and Goldberg 的 "on Rounding" 部分。
感谢我的问题收到的所有评论，我明白了发生了什么并能够解决问题。
这里我尽量总结一下。
首先，确认四舍五入不正确。如 @MarkDickinson, it can be due to a "double rounding", but I do not know if it can be confirmed. Indeed, it can also be due to others phenomenon such as the ones described in the publication given by Pascal Cuoq.
所述好像ia32 FPU在某些数字四舍五入的问题上并没有完全符合IEEE754标准。
默认情况下，GCC（32 位版本）生成使用 FPU 计算 Binary64 数字加法的代码。
但是，在我的电脑（intel core i7）上，SSE 单元也可以进行这些计算。 GCC（64 位版本）默认使用此单元。
在 GCC32 命令行上使用以下两个选项解决了我的问题。
-msse2 -mfpmath=sse.
(谢谢EOF)