几个月前写过一篇文章《log1p(x) 和 expm1(x) 函数的实现》
http://blog.csdn.net/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)的值非常小，所以这种算法的误差很小。