Goldberg的log1p与gsl_log1p
我正在寻找一个简单的log1p可移植实现。 我遇到过两个实现。
第一个出现在Theorem 4这里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) - (zx)/y ; /* cancels errors with IEEE arithmetic */ }
是否有理由偏爱另一个？
我使用log()实现测试了这两种方法，与多精度算术库相比，最大误差<0.53 ulps。 使用该log()实现作为两个log1p()变体的构建块，我发现Goldberg版本的最大错误<1.51 ulps，而GSL变体中的最大错误<0.67 ulps。 这表明后者更加准确。
在特殊情况处理方面，Goldberg变体显示出一个不匹配，因为它为+无穷大的输入返回NaN，而正确的结果是+无穷大。 对于GSL实现的特殊情况，存在三个不匹配：-1和+无穷大的输入提供NaN，而正确的结果应分别为-infinity和+ infinity。 此外，对于-0的输入，此代码返回+0，而正确的结果是-0。
如果不了解投入的分布，就很难评估绩效。 正如其他人在评论中指出的那样，当许多参数接近于零时，Goldberg的版本可能更快，因为它会跳过对这些参数的昂贵的log()调用。
两者之间没有确定的答案。 GNU科学库非常强大，使用良好，并且在所有后期版本的gcc中得到了积极支持。 使用它时，你不太可能遇到太多的惊喜。 至于你搞砸的任何其他代码，在你validation了它的逻辑并且对它的错误检查级别/方式感到满意之后，绝对没有理由不使用它。 GSL的一个小缺点是，它是您必须随身携带的另一个库，并且取决于您的代码的广泛使用将为其他平台上的其他用户提供更多挑战。 这大约是它的大小。
最好的代码是最符合项目要求的代码。