Goldberg 的 log1p 与 gsl_log1p

问题描述

我正在寻找一个简单的 log1p 可移植实现.我遇到过两种实现.

I am looking for a simple portable implementation of log1p. I have come across two implementations.

第一个在这里作为定理 4 出现http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html,

The first one appears as Theorem 4 here http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html,

上面的实现

double log1p(double p)
{
   volatile double y = p;
   return ( (1 + y) == 1 ) ? y : y * ( log( 1 + y) / ( ( 1 + y) - 1 ) );
}

第二个在 GSL http://fossies.org/dox/gsl-1.16/log1p_8c_source.html

double gsl_log1p (const double x)
{
  volatile double y, z;
  y = 1 + x;
  z = y - 1;
  return log(y) - (z-x)/y ; /* cancels errors with IEEE arithmetic */
}

是否有理由偏爱其中一种?

Is there a reason to prefer one over the other?

推荐答案

我已经使用 log() 实现测试了这两种方法，最大误差为 <0.51 ulps，与多精度算术库相比.使用 log() 实现作为两个 log1p() 变体的构建块，我发现 Goldberg 版本的最大误差为 log1p().2.5 ulps，而 GSL 变体中的最大误差小于1.5 ul.这表明后者更准确.

I have tested these two approaches using a log() implementation with a maximum error of < 0.51 ulps, comparing to a multi-precision arithmetic library. Using that log() implementation as a building block for the two log1p() variants, I found the maximum error of Goldberg's version to be < 2.5 ulps, while the maximum error in the GSL variant was < 1.5 ulps. This indicates that the latter is significantly more accurate.

在特殊情况处理方面，Goldberg 变体显示出一个不匹配，因为它为 +infinity 的输入返回 NaN，而正确的结果是 +infinity.对于 GSL 实现的特殊情况，存在三个不匹配: -1 和 +infinity 的输入传递了 NaN，而正确的结果应该分别是 -infinity 和 +infinity.此外，对于 -0 的输入，此代码返回 +0，而正确的结果是 -0.

In terms of special case handling, the Goldberg variant showed one mismatch, in that it returns a NaN for an input of +infinity, whereas the correct result is +infinity. There were three mismatches for special cases with the GSL implementation: Inputs of -1 and +infinity delivered a NaN, while the correct results should be -infinity and +infinity, respectively. Also, for an input of -0 this code returned +0, whereas the correct result is -0.

如果不了解输入的分布，就很难评估性能.正如其他人在评论中指出的那样，当许多参数接近于零时，Goldberg 的版本可能会更快，因为它跳过了对这些参数的昂贵的 log() 调用.

It is difficult to assess performance without knowledge of the distribution of the inputs. As others have pointed out in comments, Goldberg's version is potentially faster when many arguments are close to zero, as it skips the expensive call to log() for such arguments.

这篇关于Goldberg 的 log1p 与 gsl_log1p的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！