C语言 IEEE-754兼容圆形半至偶数



C标准库在C99中提供了round, lroundllround族函数。然而,这些功能不符合IEEE-754标准,因为它们不实现IEEE强制要求的"银行家四舍五入"。如果小数部分恰好为0.5,则要求将结果四舍五入到最接近的偶数值。相反,C99标准要求与零相差一半,如cppreference.com

所示。

1-3)计算最接近arg的整数值(以浮点格式),舍入到离零一半的情况,无论当前的舍入模式如何。

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

我如何实现一个函数,将四舍五入任何浮点值与半到偶数语义?如果可能的话,函数应该只依赖于语言和标准库语义,这样它就可以在非ieee浮点类型上操作。如果这是不可能的,用IEEE-754位表示定义的答案也是可以接受的。请用<limits.h><limits>描述任何常数

C标准库在C99中提供了roundlroundllround系列函数。然而,这些功能不符合IEEE-754标准,因为它们不实现IEEE强制要求的"银行家四舍五入"。

讨论单个功能是否"符合IEEE-754"是没有意义的。遵从IEEE-754要求一组具有定义语义的数据类型操作可用。它不要求这些类型或操作具有特定的名称,也不要求只有这些操作可用。实现可以提供它想要的任何附加功能,并且仍然是兼容的。如果一个实现想要提供四舍五入到奇数、四舍五入随机、四舍五入到零和不精确,它可以这样做。

IEEE-754对舍入的实际要求是提供以下六种操作:

convertToIntegerTiesToEven (x)

convertToIntegerTowardZero (x)

convertToIntegerTowardPositive (x)

convertToIntegerTowardNegative (x)

convertToIntegerTiesToAway (x)

convertToIntegerExact (x)

在C和c++中,后五个操作分别绑定到truncceilfloorroundrint函数。C11和c++ 14没有第一个绑定,但未来的版本将使用roundeven。如您所见,round实际上是必需的操作之一。

然而,roundeven在当前的实现中是不可用的,这把我们带到你问题的下一部分:

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

这个表达式的问题远远超出了"严格的IEEE-754数学"。它对负x完全不正确,对nextDown(0.5)给出错误的答案,并将2**23二进制中的所有奇数变成偶数。任何将它翻译成cvtss2si的编译器都是非常非常坏的。如果你有这样的例子,我很想看看。

我如何实现一个函数,将四舍五入任何浮点值与半到偶数语义?

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

使用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 

将数字x四舍五入,如果x与Round (x)之间的差值正好是+0.5或-0.5,并且Round (x)是奇数,则Round (x)的四舍五入方向错误,因此要从x中减去差值。

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;
}

并利用加法的自然舍入行为,而不必使用任何其他"舍入到整数"功能。

从c++ 11开始,STL中有一个函数可以进行半偶数舍入。如果将浮点舍入模式设置为FE_TONEAREST(默认值),则std::nearbyint将完成此操作。

std::nearbyint

下面是一个遵循IEEE四舍五入标准的四舍五入到偶数程序的简单实现。

逻辑:误差= 0.00001

  1. 数量= 2.5
  2. temp = floor(2.5)%2 = 2%2 = 0
  3. x = -1 + temp = -1
  4. x*error + number = 2.40009
  5. 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 for floor function)

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

给出的输出如下acc。达到标准:

(3.5)→4

(2.5)→2


编辑1:正如@Mark Dickinson在评论中指出的。可以根据需要在代码中修改错误以使其标准化。对于python,要将其转换为尽可能小的浮点值,您可以执行以下操作。

import sys
error = sys.float_info.min

最新更新