log1p(x) 函数的实现 (续)

几个月前写过一篇文章《log1p(x) 和 expm1(x) 函数的实现》

http://blog.youkuaiyun.com/liyuanbhu/article/details/8544644

当时给出了个 GSL 的 log1p(x) 的代码实现。不过当时没想明白为什么要那样算，最近在水木上看一个帖子时突然就领悟了。

GSL 上的代码如下：

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 */  
}  

我们知道计算机上的浮点数是只有有限的精度的。

所以 y:=1+x 这条赋值语句当x较小时算出的y并不精确等于1+x。

或者可以这样认为，我们实际计算的是

y:1+x'

而x’ 就是计算1+x 时小数点对其的过程种x有效位数损失后的结果。 

因此，直接计算 log(1+x) 实际上是计算了 log(1+x')。

由于 x 与 x' 很接近，所以下面的式子近似成立：

上面的式子变一下型就得到：

这也就是GSL 上采用的代码了(x’ 对应代码中的z)。

多说一句，在C程序中不能写为：

z= 1+x-1;

否则编译器会自作聪明的优化成 z = x，并且似乎无法通过编译器的编译选项关掉这种优化。

另外，有一篇很有名的文章，What Every Computer Scientist Should Know About Floating-Point Arithmetic，上面也探讨了这个计算。并且给出了另一个计算方法。

这个算法的依据在于当 x 与 x’ 很接近时，

x(x’-x)的值非常小，所以这种算法的误差很小。