log1p(x) 函數的實現 (續)

當時給出了個 GSL 的 log1p(x) 的代碼實現。不過當時沒想明白為什麼要那樣算，最近在水木上看一個帖子時突然就領悟了。

GSL 上的代碼如下：

[cpp] view plaincopyprint?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 */   
}   

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)的值非常小，所以這種算法的誤差很小。