3

问题

对于实现精确 IEEE 754 算术的 C99 编译器f，divisor类型 的值是否float存在使得f / divisor != (float)(f * (1.0 / divisor))?

编辑：通过“实现精确的 IEEE 754 算术”，我的意思是一个正确地将 FLT_EVAL_METHOD 定义为 0 的编译器。

语境

提供符合 IEEE 754 的浮点的 AC 编译器只能将单精度除以常数替换为单精度乘以逆，如果所述逆本身可以精确表示为float.

在实践中，这只发生在二次幂的情况下。因此，程序员 Alex 可能确信它f / 2.0f会像以前一样被编译f * 0.5f，但如果 Alex 可以接受乘以0.10f而不是除以 10，则 Alex 应该通过在程序中编写乘法来表达它，或者使用编译器选项，例如 GCC 的-ffast-math.

这个问题是关于将单精度除法转换为双精度乘法。它总是产生正确的舍入结果吗？它是否有可能更便宜，从而成为编译器可能进行的优化（即使没有-ffast-math）？

我已经比较了 1 到 2(float)(f * 0.10)之间f / 10.0f的所有单精度值f，但没有找到任何反例。float这应该涵盖产生正常结果的 normal 的所有划分。

然后我使用以下程序将测试推广到所有除数：

#include <float.h>
#include <math.h>
#include <stdio.h>

int main(void){
  for (float divisor = 1.0; divisor != 2.0; divisor = nextafterf(divisor, 2.0))
    {
      double factor = 1.0 / divisor; // double-precision inverse
      for (float f = 1.0; f != 2.0; f = nextafterf(f, 2.0))
        {
          float cr = f / divisor;
          float opt = f * factor; // double-precision multiplication
          if (cr != opt)
            printf("For divisor=%a, f=%a, f/divisor=%a but (float)(f*factor)=%a\n",
                   divisor, f, cr, opt);
        }
    }
}

搜索空间足够大，足以让这变得有趣 (2 46 )。该程序当前正在运行。有人可以告诉我它是否会在完成之前打印一些东西，也许会解释为什么或为什么不打印？

4

3 回答 3

6

你的程序不会打印任何东西，假设是圆角到偶数的舍入模式。论证的实质如下：

我们假设两者f和divisor都在1.0和之间2.0。所以对于一些整数f = a / 2^23和在范围内。这个案例并不有趣，所以我们可以进一步假设。divisor = b / 2^23ab[2^23, 2^24)divisor = 1.0b > 2^23

(float)(f * (1.0 / divisor))可能给出错误结果的唯一方法是使精确值f / divisor非常接近中间情况（即，两个单精度浮点数之间正好中间的数字），以至于表达式中的累积错误f * (1.0 / divisor)将我们推到另一边从真正的价值中途的情况。

但这不可能发生。为简单起见，我们首先假设f >= divisor, 使得精确商在[1.0, 2.0). 现在，间隔中单精度的任何中途情况都有一些奇整数的[1.0, 2.0)形式。的确切值是，因此差的绝对值在下界，因此至少是（因为）。因此，我们距离任何中途情况都不止双精度 ulps，而且应该很容易证明双精度计算中的误差永远不会超过 16 ulps。（我没有做算术，但我猜很容易在错误上显示 3 ulps 的上限。）c / 2^24c2^24 < c < 2^25f / divisora / bf / divisor - c / 2^241 / (2^24 b)1 / 2^48b < 2^2416

所以f / divisor不能足够接近中途案例来制造问题。请注意，这f / divisor也不能是一个精确的中途情况：因为c是奇数，c并且2^24是互质的，所以我们唯一可以拥有的方法c / 2^24 = a / b是 ifb是 的倍数2^24。但是b在范围内(2^23, 2^24)，所以这是不可能的。

类似的情况f < divisor：中间情况具有形式c / 2^25，并且类似的论点表明abs(f / divisor - c / 2^25)大于1 / 2^49，这再次为我们提供了16双精度 ulps 的余量。

于 2013-10-25T20:16:33.140 回答
2

如果非默认舍入模式是可能的，那肯定是不可能的。例如，在替换3.0f / 3.0f为 时3.0f * C，C小于精确倒数的值将在向下或接近零的舍入模式下产生错误的结果，而C大于精确倒数的值将在向上舍入模式下产生错误的结果。

如果您限制为默认舍入模式，我不太清楚您正在寻找的内容是否可能。如果我有任何想法，我会考虑并修改这个答案。

于 2013-10-25T12:55:13.623 回答
1

随机搜索产生了一个例子。

看起来当结果是“非正规/次正规”数字时，不等式是可能的。但是，也许我的平台不符合 IEEE 754 标准？

f        0x1.7cbff8p-25 
divisor -0x1.839p+116
q       -0x1.f8p-142
q2      -0x1.f6p-142

int MyIsFinite(float f) {
  union {
    float f;
    unsigned char uc[sizeof (float)];
    unsigned long ul;
  } x;
  x.f = f;
  return (x.ul & 0x7F800000L) != 0x7F800000L;
}

float floatRandom() {
  union {
    float f;
    unsigned char uc[sizeof (float)];
  } x;
  do {
    size_t i;
    for (i=0; i<sizeof(x.uc); i++) x.uc[i] = rand();
  } while (!MyIsFinite(x.f));
  return x.f;
}

void testPC() {
  for (;;) {
    volatile float f, divisor, q, qd;
    do {
      f = floatRandom();
      divisor = floatRandom();
      q = f / divisor;
    } while (!MyIsFinite(q));
    qd = (float) (f * (1.0 / divisor));
    if (qd != q) {
      printf("%a %a %a %a\n", f, divisor, q, qd);
      return;
    }
  }

}

Eclipse PC 版本：Juno Service Release 2 Build id：20130225-0426

于 2013-10-25T14:49:41.007 回答