符合 IEEE-754 的舍入到偶数答案

【问题标题】：IEEE-754 compliant round-half-to-even符合 IEEE-754 的舍入到偶数
【发布时间】：2015-12-21 04:21:15
【问题描述】：

C 标准库在 C99 中提供了 round、lround 和 llround 系列函数。但是，这些函数不符合 IEEE-754 标准，因为它们没有按照 IEEE 的要求实现半对偶的“银行家四舍五入”。如果小数部分正好是 0.5，则半到偶数舍入要求将结果舍入到最接近的偶数值。正如cppreference.com 中所述，C99 标准要求从零开始一半

1-3) 计算最接近 arg 的整数值（以浮点格式），从零开始舍入一半的情况，而不管当前的舍入模式。

在 C 中实现舍入的常用方式是表达式 (int)(x + 0.5f)，尽管它在严格的 IEEE-754 数学中是 incorrect，但通常由编译器翻译成正确的 cvtss2si 指令。但是，这当然不是一个可移植的假设。

如何实现一个函数，该函数将使用半偶数语义舍入任何浮点值？如果可能，该函数应该只依赖于语言和标准库语义，以便它可以对非 IEEE 浮点类型进行操作。如果这是不可能的，那么根据 IEEE-754 位表示定义的答案也是可以接受的。请用<limits.h> 或<limits> 来描述任何常量。

【问题讨论】：

它们符合标准，但没有全部实现。 “四舍五入，远离零”是 IEEE-754 舍入模式之一。
IEEE-754 指出存在教科书舍入模式，但仅适用于十进制浮点数。 An implementation of this standard shall provide roundTiesToEven and the three directed rounding attributes. A decimal format implementation of this standard shall provide roundTiesToAway as a userselectable rounding-direction attribute. The rounding attribute roundTiesToAway is not required for a binary format implementation. 将关系四舍五入到偶数是唯一可移植（且合理）的模式。
而不是round()，请查看rint()，默认舍入模式为“舍入到最接近或偶数”。
您可能对 CPython 的解决方案感兴趣，here。
@njuffa 非常感谢指向rint（非常不直观的名称）的指针。我不得不仔细检查 C 标准，但是当设置 __STDC_IEC_559__ 时，它似乎确实做了正确的事情。真正的谜团是为什么round 做的事情与rint 和FE_TONEAREST 不同。

标签： c++ c floating-point ieee-754 bankers-rounding

【解决方案1】：

将数字 x 舍入，如果 x 和舍入 (x) 之间的差值正好是 +0.5 或 -0.5，并且舍入 (x) 是奇数，那么舍入 (x) 的舍入方向错误，所以你减去与 x 的区别。

【讨论】：

如果 x 的小数部分正好是一半，是否保证 (round(x) - x) 将返回 0.5/-0.5，即这是一个精确的操作？
是的。 round(x) - x 的结果将是完全可表示的（假设为 IEEE 754 格式）。
一个邪恶的技巧是在fabs(x - round(x)) == 0.5 的情况下使用2.0 * round(0.5 * x)，这样就省去了判断round(x) 是偶数还是奇数的有点混乱的任务。

【解决方案2】：

使用 C 标准库中的 remainder(double x, 1.0)。这与当前的舍入模式无关。

余数函数计算IEC 60559所需的余数x REM y

remainder() 在这里很有用，因为它满足了 OP 与偶数要求的联系。

double round_to_nearest_ties_to_even(double x) {
  x -= remainder(x, 1.0);
  return x;
}

测试代码

void rtest(double x) {
  double round_half_to_even = round_to_nearest_ties_to_even(x);
  printf("x:%25.17le   z:%25.17le \n", x, round_half_to_even);
}

void rtest3(double x) {
  rtest(nextafter(x, -1.0/0.0));
  rtest(x);
  rtest(nextafter(x, +1.0/0.0));
}

int main(void) {
  rtest3(-DBL_MAX);
  rtest3(-2.0);
  rtest3(-1.5);
  rtest3(-1.0);
  rtest3(-0.5);
  rtest(nextafter(-0.0, -DBL_MAX));
  rtest(-0.0);
  rtest(0.0);
  rtest(nextafter(0.0, +DBL_MAX));
  rtest3(0.5);
  rtest3(1.0);
  rtest3(1.5);
  rtest3(2.0);
  rtest3(DBL_MAX);
  rtest3(0.0/0.0);
  return 0;
}

输出

x:                     -inf   z:                     -inf 
x:-1.79769313486231571e+308   z:-1.79769313486231571e+308 
x:-1.79769313486231551e+308   z:-1.79769313486231551e+308 
x: -2.00000000000000044e+00   z: -2.00000000000000000e+00 
x: -2.00000000000000000e+00   z: -2.00000000000000000e+00 
x: -1.99999999999999978e+00   z: -2.00000000000000000e+00 
x: -1.50000000000000022e+00   z: -2.00000000000000000e+00 
x: -1.50000000000000000e+00   z: -2.00000000000000000e+00 tie to even
x: -1.49999999999999978e+00   z: -1.00000000000000000e+00 
x: -1.00000000000000022e+00   z: -1.00000000000000000e+00 
x: -1.00000000000000000e+00   z: -1.00000000000000000e+00 
x: -9.99999999999999889e-01   z: -1.00000000000000000e+00 
x: -5.00000000000000111e-01   z: -1.00000000000000000e+00 
x: -5.00000000000000000e-01   z:  0.00000000000000000e+00 tie to even 
x: -4.99999999999999944e-01   z:  0.00000000000000000e+00 
x:-4.94065645841246544e-324   z:  0.00000000000000000e+00 
x: -0.00000000000000000e+00   z:  0.00000000000000000e+00 
x:  0.00000000000000000e+00   z:  0.00000000000000000e+00 
x: 4.94065645841246544e-324   z:  0.00000000000000000e+00 
x:  4.99999999999999944e-01   z:  0.00000000000000000e+00 
x:  5.00000000000000000e-01   z:  0.00000000000000000e+00 tie to even 
x:  5.00000000000000111e-01   z:  1.00000000000000000e+00 
x:  9.99999999999999889e-01   z:  1.00000000000000000e+00 
x:  1.00000000000000000e+00   z:  1.00000000000000000e+00 
x:  1.00000000000000022e+00   z:  1.00000000000000000e+00 
x:  1.49999999999999978e+00   z:  1.00000000000000000e+00 
x:  1.50000000000000000e+00   z:  2.00000000000000000e+00 tie to even 
x:  1.50000000000000022e+00   z:  2.00000000000000000e+00 
x:  1.99999999999999978e+00   z:  2.00000000000000000e+00 
x:  2.00000000000000000e+00   z:  2.00000000000000000e+00 
x:  2.00000000000000044e+00   z:  2.00000000000000000e+00 
x: 1.79769313486231551e+308   z: 1.79769313486231551e+308 
x: 1.79769313486231571e+308   z: 1.79769313486231571e+308 
x:                      inf   z:                      inf 
x:                      nan   z:                      nan 
x:                      nan   z:                      nan 
x:                      nan   z:                      nan

【讨论】：

【解决方案3】：

C 标准库在 C99 中提供了 round、lround 和 llround 系列函数。但是，这些函数不符合 IEEE-754 标准，因为它们没有按照 IEEE 的要求实现对半偶数的“银行家四舍五入”...

谈论单个函数是否“符合 IEEE-754”是没有意义的。 IEEE-754 合规性要求一组具有定义语义的数据类型操作可用。它不要求这些类型或操作具有特定名称，也不要求仅这些操作可用。一个实现可以提供它想要的任何附加功能并且仍然是合规的。如果实现想要提供奇数舍入、随机舍入、零舍入和不精确时陷阱，则可以这样做。

IEEE-754 对四舍五入的实际要求是提供了以下六个操作：

convertToIntegerTiesToEven(x)

convertToIntegerTowardZero(x)

convertToIntegerTowardPositive(x)

convertToIntegerTowardNegative(x)

convertToIntegerTiesToAway(x)

convertToIntegerExact(x)

在 C 和 C++ 中，最后五个操作分别绑定到 trunc、ceil、floor、round 和 rint 函数。 C11 和 C++14 第一个没有绑定，但未来的修订版将使用 roundeven。如您所见，round 实际上是必需的操作之一。

但是，roundeven 在当前实现中不可用，这将我们带到您问题的下一部分：

在 C 中实现舍入的常用方式是表达式 (int)(x + 0.5f)，尽管它在严格的 IEEE-754 数学中不正确，但通常由编译器翻译成正确的 cvtss2si 指令。但是，这当然不是一个可移植的假设。

该表达式的问题远远超出了“严格的 IEEE-754 数学”。否定的x 是完全不正确的，给出了nextDown(0.5) 的错误答案，并将 2**23 binade 中的所有奇数转换为偶数。任何将它翻译成cvtss2si 的编译器都被严重破坏了。如果你有这样的例子，我很乐意看到。

如何实现一个函数，该函数将使用半偶数语义舍入任何浮点值？

正如 njuffa 在评论中指出的那样，您可以确保设置了默认舍入模式并使用rint（或lrint，因为听起来您实际上想要一个整数结果），或者您可以通过调用round 来实现自己的舍入函数，然后像 gnasher729 建议的那样修复中途情况。一旦采用了 C 的 n1778 绑定，您就可以使用 roundeven 或 fromfp 函数来执行此操作，而无需控制舍入模式。

【讨论】：

我同意 (int)(x + 0.5f) 是错误的，原因有很多。但是，MSVC 特别检测到这个习语并用 cvtss2si 替换它，我猜这就是它常用的原因。
@68ejxfcj5669：我怀疑是这样，但我很想被证明是错误的。我没有要测试的 Windows 机器，在网上搜索也没有发现任何关于这种转变的说法。您有 MSVC 执行此操作的程序的实际示例吗？

【解决方案4】：

float 数据类型可以表示 8388608.0f 到 16777216.0f 范围内的所有整数，但不能表示分数。任何大于 8388607.5f 的 float 数字都是整数，不需要四舍五入。将 8388608.0f 添加到任何小于该值的非负 float 将产生一个整数，该整数将根据当前的舍入模式（通常是四舍五入到偶数）进行四舍五入。然后减去 8388608.0f 将产生一个正确舍入的原始版本（假设它在合适的范围内）。

因此，应该可以执行以下操作：

float round(float f)
{
  if (!(f > -8388608.0f && f < 8388608.0f)) // Return true for NaN
    return f;
  else if (f > 0)
    return (float)(f+8388608.0f)-8388608.0f;
  else
    return (float)(f-8388608.0f)+8388608.0f;
}

并利用加法的自然舍入行为，而无需使用任何其他“舍入到整数”工具。

【讨论】：

【解决方案5】：

下面是round-half to even程序的简单实现，它遵循IEEE的舍入标准。

逻辑： 误差 = 0.00001

数字=2.5

temp = floor(2.5)%2 = 2%2 = 0

x = -1 + temp = -1

x*error + number = 2.40009

round(2.40009) = 2

注意：这里的误差是0.00001，即如果出现2.500001，那么它会四舍五入到2而不是3

Python 2.7 实现：

temp = (number)     
rounded_number = int( round(-1+ temp%2)*0.00001 + temp )

C++ 实现：（使用 math.h 实现 floor 函数）

float temp = (number)     
int rounded_number = (int)( (-1+ temp%2)*0.00001 + temp + 0.5)

这将给出的输出如下所示。符合标准：

(3.5) -> 4

(2.5) -> 2

编辑 1： 正如@Mark Dickinson 在 cmets 中指出的那样。可以根据代码中的要求修改错误以使其标准化。对于python，要将其变成可能的最小浮点值，您可以执行以下操作。

import sys
error = sys.float_info.min

【讨论】：

如果 2.500001 舍入为 2 而不是 3，那么这不是符合标准的舍入。
该程序是为希望在自己的代码中实现此功能的人编写的。错误可以更改为代码中所需的任何精度点，即，它甚至可以是可能的最小浮点值。 error = sys.float_info.min（对于 Python），以及类似的 C++ 等价物
好的，但它没有回答问题：OP 正在寻找符合 IEEE 754 的舍入操作。这个答案中的代码没有给出。并且使用error = sys.float_info.min 不是解决方案。如果您测试您的代码，您会发现即使是当前代码也不会从3.5 的输入中给出4；它给了3。

【解决方案6】：

自 C++11 起，STL 中有一个函数可以进行半偶数舍入。如果浮点舍入模式设置为FE_TONEAREST（默认），那么std::nearbyint 就可以了。

C++ Reference for std::nearbyint

【讨论】：