实施单precision师双链precision乘法 [英] Implementing single-precision division as double-precision multiplication

问题描述

对于C99编译器实现精确的IEEE 754算法，做的f的数值，除数类型浮动存在这样 F /除数！=（浮点）（F *（1.0 /除数））？

For a C99 compiler implementing exact IEEE 754 arithmetic, do values of f, divisor of type float exist such that f / divisor != (float)(f * (1.0 / divisor))?

编辑：由实施精确的IEEE 754算法我的意思是编译器，它理所当然地FLT_EVAL_METHOD定义为0

By "implementing exact IEEE 754 arithmetic" I mean a compiler that rightfully defines FLT_EVAL_METHOD as 0.

交流编译如果所述逆本身重新presentable提供IEEE 754标准的浮点只能由一个常数由逆的单precision乘法取代的单precision除法正是因为一个浮动。

A C compiler that provides IEEE 754-compliant floating-point can only replace a single-precision division by a constant by a single-precision multiplication by the inverse if said inverse is itself representable exactly as a float.

在实践中，这仅发生于两个权力。因此，一个程序员，亚历克斯，可以确信 F / 2.0F 将被编译，就好像它一直 F * 0.5F ，但如果它是可以接受的亚历克斯通过 0.10f 繁衍，而不是10分，亚历克斯应该在程序编写乘法前preSS它，或通过使用编译器选项，如GCC的 -ffast-数学。

In practice, this only happens for powers of two. So a programmer, Alex, may be confident that f / 2.0f will be compiled as if it had been f * 0.5f, but if it is acceptable for Alex to multiply by 0.10f instead of dividing by 10, Alex should express it by writing the multiplication in the program, or by using a compiler option such as GCC's -ffast-math.

这个问题是关于转化单precision划分为一个双precision倍增。它总是产生正确舍入的结果？是否有一个机会，它可能是更便宜，因而是一种优化编制人员可能使（即使没有 -ffast-数学）？

This question is about transforming a single-precision division into a double-precision multiplication. Does it always produce the correctly rounded result? Is there a chance that it could be cheaper, and thus be an optimization that compilers might make (even without -ffast-math)?

我比较（浮点）（F * 0.10）和 F / 10.0f 所有单$ ˚F 1和2之间的的对$ pcision值都没找到反例。这应包括正常浮动所有部门取值产生一个正常的结果。

I have compared (float)(f * 0.10) and f / 10.0f for all single-precision values of f between 1 and 2, without finding any counter-example. This should cover all divisions of normal floats producing a normal result.

然后，我概括了测试与下面的程序中的所有除数：

Then I generalized the test to all divisors with the program below:

#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 ⁴⁶ ）。该方案目前正在运行。谁能告诉我它是否会打印一些东西，也许解释为什么或者为什么不呢，它已经完成过吗？

The search space is just large enough to make this interesting (2⁴⁶). The program is currently running. Can someone tell me whether it will print something, perhaps with an explanation why or why not, before it has finished?

推荐答案

您的程序不会打印任何东西，假设圆关系到偶数舍入模式。争论的本质是如下：

Your program won't print anything, assuming round-ties-to-even rounding mode. The essence of the argument is as follows:

我们假设这两个˚F和除数之间是 1.0 和 2.0 。因此， F = A / 2 ^ 23 和除数= B / 2 ^ 23 对于一些整数在和 b 范围 [2 ^ 23，2 ^ 24）。此案除数= 1.0 不感兴趣，所以我们可以进一步假设 B＆GT; 2 ^ 23 。

We're assuming that both f and divisor are between 1.0 and 2.0. So f = a / 2^23 and divisor = b / 2^23 for some integers a and b in the range [2^23, 2^24). The case divisor = 1.0 isn't interesting, so we can further assume that b > 2^23.

的唯一方法是（浮点）（F *（1.0 /除数））可以给错误的结果将是精确值˚F /除数是如此接近一半的情况下（即数正中间两个单核precision浮筒之间），在恩pression的累积误差 F *（1.0 /除数）我们推到与真实值一半案件的对方的

The only way that (float)(f * (1.0 / divisor)) could give the wrong result would be for the exact value f / divisor to be so close to a halfway case (i.e., a number exactly halfway between two single-precision floats) that the accumulated errors in the expression f * (1.0 / divisor) push us to the other side of that halfway case from the true value.

但是，这是不可能发生的。为简单起见，我们首先假设 F＆GT; =除数，以便准确的商在 [1.0,2.0）。现在在区间单precision任何中途情况下 [1.0,2.0）的形式 C / 2 ^ 24 一些奇怪的整数 C 与 2 ^ 24 LT; C＆LT; 2 ^ 25 。 的F /除数的精确值 A / B ，所以差的绝对值 F /除数 - C / 2 ^ 24 是 1 /（2 ^ 24 b），所以至少 1 / ^ 48 （因为 b 2 ^ 24 ）。所以我们超过 16 双击precision ULPS任何中途离开的情况下，它应该很容易显示，在双precision错误计算不能超过16 ULPS。 （我没有做算术，但我猜这是容易证明的上限对错误3 ULPS。）

But that can't happen. For simplicity, let's first assume that f >= divisor, so that the exact quotient is in [1.0, 2.0). Now any halfway case for single precision in the interval [1.0, 2.0) has the form c / 2^24 for some odd integer c with 2^24 < c < 2^25. The exact value of f / divisor is a / b, so the absolute value of the difference f / divisor - c / 2^24 is bounded below by 1 / (2^24 b), so is at least 1 / 2^48 (since b < 2^24). So we're more than 16 double-precision ulps away from any halfway case, and it should be easy to show that the error in the double precision computation can never exceed 16 ulps. (I haven't done the arithmetic, but I'd guess it's easy to show an upper bound of 3 ulps on the error.)

所以 F /除数不能足够接近一半的情况下创造的问题。需要注意的是 F /除数不能是的确切的一半的情况下，无论是：因为 C 是奇数， C 和 2 ^ 24 互质，所以唯一的办法，我们可以有 C / 2 ^ 24 = A / b 是，如果 b 是一个多 2 ^ 24 。但 B 的范围为（2 ^ 23，2 ^ 24），所以这是不可能的。

So f / divisor can't be close enough to a halfway case to create problems. Note that f / divisor can't be an exact halfway case, either: since c is odd, c and 2^24 are relatively prime, so the only way we could have c / 2^24 = a / b is if b is a multiple of 2^24. But b is in the range (2^23, 2^24), so that's not possible.

的情况下 F＆LT;除数是相似的：在半路的情况下则有如下形式 C / 2 ^ 25 和类似的说法表明， ABS（ F /除数 - C / 2 ^ 25）大于 1 / ^ 49 ，这再次给了我们一个保证金 16 双击precision ULPS一起玩。

The case where f < divisor is similar: the halfway cases then have the form c / 2^25 and the analogous argument shows that abs(f / divisor - c / 2^25) is greater than 1 / 2^49, which again gives us a margin of 16 double-precision ulps to play with.

这篇关于实施单precision师双链precision乘法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！