在STDC中究竟做了什么(1.0e300 + pow(2.0,-30.0)> 1.0)?

问题描述 投票:8回答:1

我遇到了一个计算atan(x)的函数(源代码是here)。将它减少到我的问题的核心并稍微重新格式化,他们有类似的东西:

static const double one   = 1.0,
                   huge   = 1.0e300;

double atan(double x)
{
  /* A lot of uninteresting stuff here */

  if (ix < 0x3fdc0000) {              /* |x| < 0.4375 */

    if (ix < 0x3e200000) {            /* |x| < 2^-29 */
      if ((huge + x) > one) return x; /* raise inexact */
    }

    id = -1;
  }

  /* A lot of more uninteresting stuff */
}

我会对if ((huge + x) ...应该做什么以及如何工作非常感兴趣。

根据评论,如果x的绝对值小于2^-29,表达式或比较会引发inexact错误。

我的第一个问题是,我目前无法理解为什么这样做:如果arctan的绝对值太小,使用该函数计算x会导致不精确的结果,为什么他们不使用像if (fabs(x) < [some_value_here]) ...这样的东西?我怀疑这只是因为inexact警告不会在他们的硬件/库中以这种方式提出,但我想知道肯定。

假设我是对的,我的第二个问题是我不明白为什么需要比较。我认为这里的关键点是非常小的数字被添加到一个非常大的数字中,因此这种添加不会充分改变大数字,甚至根本不会改变。因此,它是增加inexact警告,而不是比较。所以我问自己比较应该做些什么。这只是为了强制编译器实际计算(huge + x),否则可能会进行优化?

最后,如果有人能够解释一下数学,我将不胜感激。为1.0e300选择值huge似乎是一个相当随意的选择。但这只是一个额外的问题,因为我承认我还没有完成我的作业的数学部分(我不是关于double值及其IEEE754表示的全新手,但是理解这段代码的数学方面需要一些时间,除非有人给出简短的解释)。

编辑1

偶然看到:

该函数的float32版本,包括上面讨论的怪异线,几乎完全仍然在glibc 2.19中!由于glibc是可移植的,因此该代码也应如此。它位于子目录sysdeps\ieee754\flt-32中,所以我想这是float32函数的软件仿真,其中可移植性没有问题,因为硬件相关的奇怪性不会出现(我认为软件仿真完全按照IEEE754)。

c ieee-754 atan2
1个回答
7
投票

if ((huge + x) > one) return x;的目的是生成浮点不精确异常,然后从例程返回。

浮点异常不是陷阱或处理器异常。它只是意味着在浮点运算中发生了一些不寻常的事情。然后会发生什么取决于操作的情况。特别是,可以设置浮点环境,以便不精确的异常仅在特殊寄存器中引发标志并继续操作,从而提供数值结果。或者可以设置它以便不精确的异常导致陷阱,并且程序控制被重定向到陷阱处理程序。

实现atan的代码不知道如何设置浮点环境。也许它可以获取设置,但它不想打扰它。鉴于它已经确定无法精确计算反正切函数,它触发浮点不精确异常的最简单方法仅仅是执行一个具有不精确结果的简单加法。这种不精确的添加将具有与不精确的反正切所需的相同行为 - 它将仅提升标志或导致陷阱,具体取决于设置。

至于为何与ix < 0x3e200000进行比较,目前尚不清楚。一方面,ix已被调整以反映绝对值,而x没有,所以为什么不使用已经准备好的ix而不是使用另一种操作来生产fabs(x)?此外,整数比较通常比浮点比较占用更少的处理器资源,尤其是在编写此代码时的处理器中。或者可能是作者只是在作者身上使用了一个,也许使用ix编写大部分代码来操作浮点编码而不是x来操作浮点值,他们不想切换来回不必要的。也可能是因为代码是在十六进制浮点表示法可用之前编写的(所以我们可以编写x < 0x1p-29f),并且编译器不能很好地将十进制数字转换为浮点值,所以他们不想写出来源代码中的浮点值。

这种代码是有问题的,并且高度依赖于它所编写的C实现。通常,编译器可能无法保证在程序执行期间会评估(huge + x) > one。编译器可能会在编译时对其进行评估。但是,据推测,如果此代码是为特定的C实现编写的,则他们知道编译器将在编译时对其进行评估,或者确保实现相同的结果,包括引发浮点不精确异常。

从表面上看,(huge + x) > one似乎没有做任何有用的huge + x单独做不到,但也许作者知道我们没有的编译器。

huge不需要是1.0e300。任何这么大的值,以及hugex的总和都不够精确就足够了。

© www.soinside.com 2019 - 2024. All rights reserved.