将单精度除法实现为双精度乘法

Question

对于实现精确 IEEE 754 算术的 C99 编译器，是否存在浮点型的 f、除数的值，使得 f / 除数 ！= （浮点数）（f * （1.0 / 除数））） ？

编辑：所谓“实现精确的IEEE754算法”，我指的是一个正确地将FLT_EVAL_METHOD定义为0的编译器。

提供符合IEEE 754标准的浮点数的C编译器只能将单精度除法替换为常数，如果所述逆本身可以完全表示为浮点数。

实际上，这种情况只发生在2的幂上。因此，程序员Alex可能确信< code>f / 2.0f将被编译为< code>f * 0.5f，但是如果Alex可以乘以< code>0.10f而不是除以10，那么Alex应该通过在程序中编写乘法来表示，或者通过使用诸如GCC的< code>-ffast-math这样的编译器选项来表示。

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

我比较了（浮点数）（f*0.10）和f/10.0f在1和2之间的f的所有单精度值，没有找到任何反例。这应该涵盖产生正常结果的正常浮点数的所有划分。

然后我用下面的程序将测试推广到所有的除数:

#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⁴⁶）。程序当前正在运行。在完成之前，有人能告诉我它是否会打印一些东西，也许会解释为什么或为什么不打印？

Answer 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内部版本id:20130225-0426

Answer 2

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

我不太清楚如果你限制在默认的舍入模式下，你所寻找的是否可能。如果我想出什么，我会考虑一下并修改这个答案。

Answer 3

您的程序不会打印任何内容，假设是舍入-绑扎到偶数舍入模式。该论点的本质如下：

我们假设 f 和除数都在 1.0 和 2.0 之间。因此，对于 [2^23， 2^24） 范围内的某些整数 a 和 b，f = a / 2^23 且除数 = b / 2^23。大小写除数 = 1.0 并不有趣，因此我们可以进一步假设 b

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

但这不可能发生。为了简单起见，让我们首先假设< code>f

因此，f /除数不能足够接近中途情况来产生问题。请注意，f /除数也不可能是精确的中途情况：由于c是奇数，c和2 ^ 24是相对质数的，所以我们唯一可以有c / 2 ^ 24 = a / b的方法是如果b是2 ^ 24的倍数。但是 b 在 （2^23， 2^24） 的范围内，所以这是不可能的。

f的情况