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

Question

题

为了实现精确的IEEE 754算术C99编译器,这样做的价值f,divisor类型的float存在,使得f / divisor != (float)(f * (1.0 / divisor))？

编辑:通过"实现精确的IEEE 754算术"我的意思是一个正确定义FLT_EVAL_METHOD为0的编译器.

上下文

提供符合IEEE 754标准的浮点的AC编译器只能通过单次精度乘以逆来用常量替换单精度除法,如果所述逆本身可以完全表示为a float.

实际上,这只发生在两个人的权力上.因此,一个程序员Alex可能会确信它f / 2.0f会被编译为好像f * 0.5f,但是如果Alex可以接受0.10f而不是除以10,那么Alex应该通过在程序中编写乘法来表达它,或者使用一个编译器选项,如GCC -ffast-math.

这个问题是关于将单精度除法转换为双精度乘法.它总能产生正确的圆形结果吗？它是否有可能更便宜,因此是编译器可能做出的优化(即使没有-ffast-math)？

我比较(float)(f * 0.10)和f / 10.0f为所有单精度值f1和2之间,没有发现任何反例.这应该涵盖float正常结果的所有分区,产生正常结果.

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

#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

假设圆形连接到均匀舍入模式,您的程序将不会打印任何内容.论证的实质如下:

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

唯一(float)(f * (1.0 / divisor))可能产生错误结果的方法是使精确值f / divisor接近中间情况(即两个单精度浮点数之间的数字),表达式中的累积误差f * (1.0 / divisor)将我们推向另一侧那个中间案例的真实价值.

但这不可能发生.为简单起见,我们首先假设,f >= divisor确切的商是这样的[1.0, 2.0).现在对于在区间单精度任何中途的情况下[1.0, 2.0)具有如下形式c / 2^24对于一些奇整数c用2^24 < c < 2^25.确切的值f / divisor是a / b,所以差值的绝对值f / divisor - c / 2^24低于下限1 / (2^24 b),因此至少是1 / 2^48(因为b < 2^24).因此,我们不得超过16任何中途情况下的双精度ulps,并且应该很容易证明双精度计算中的误差永远不会超过16 ulps.(我没有做过算术,但是我猜想在错误上显示3 ulps的上限很容易.)

所以f / divisor不能足够接近中途案例来制造问题.注意,f / divisor不能是一个确切的中间情况,因为:因为c是奇数,c并且2^24是相对素数,所以我们唯一可能的方法c / 2^24 = a / b是if b是倍数2^24.但是b在范围内(2^23, 2^24),所以这是不可能的.

类似的情况f < divisor:中途案例然后有形式c / 2^25,类似的论证表明abs(f / divisor - c / 2^25)大于1 / 2^49,这再次给我们一个16双精度ulps 的余地.