[C語言數值分析] cmath / math.h - log (nature log)－Edison.X. Blog

nature log wiki 網頁，相關的函式還有 log10。

求 nature log 可以用積分法，如 ln(x) = [ 1/t ]dt ， t E [1, x]，

但這並不實際，速度反而慢。簡單的方式是用下列級數。

ln(x) = (x-1) - (x-1)^2/2 + (x-1)^3/3 -..., for 0<x<2

但這級數收斂速度頗慢，且 x 有範圍限制，並不是很好用的級數，有興趣可先實做。

而另一組公式是先令 y = (x-1) / (x+1)，

ln(x) = ln ( (1+y) / (1-y) ) = 2y * ( 1 + y^2 / 3 + y^4 / 5 +... + y^2n / (2n+1) )， for x > 0

基於此，程式碼如下示 ( eps 已於函式外做過定義 )

double xln_small(const double x)
{
    const double m=(x-1.0)/(x+1.0);
    const double m2 = m*m;
    double ans1, ans2=m, pwr=m, cnt=1.0;
    do{
        ans1=ans2;
        pwr*=m2, cnt+=2.0;
        ans2+=(pwr/cnt);
    }while(fabs(ans1-ans2)>eps);
    return 2.0 * ans2;
}

注意了，這份 code 我有加上 _small 修飾，

原因是在 x 值大的時候，收斂速度整個慢下來。

為了加快收斂之速度，用了一點 nature log 之特性。

假設輸入引數為 x ，將其分解成 x = A * 2^B，

ln ( x ) = ln ( A * 2^B ) = ln (A) + B * ln (2)。

上面看起來是要求二個 nature log，實際上 ln( 2 ) 可以直接先存起來，

到時候直接調用即可。至於 x = A * 2^B 該怎麼做？答案是：去分析 IEEE754。

這裡偷懶，直接調用 math.h 裡之 frexp 完成。程式碼參考如下。

double xln(const double x) // slow, 355 / 1219
{    
    int B;
    const double A = frexp(x, &B); // x = A * 2^B

    // find ln(A) = 2.0*ans2
    const double m=(A-1.0)/(A+1.0);
    const double m2 = m*m;
    double ans1, ans2=m, pwr=m, cnt=1.0;
    do{
        ans1=ans2;
        pwr*=m2, cnt+=2.0;
        ans2+=(pwr/cnt);
    }while(fabs(ans1-ans2)>eps);
    // return ln(A) + B*ln(2)
    /* LN_2 has defined as 6.9314718055994529E-1*/
    return 2.0*ans2 + B*LN_2; 
}

和單單使用公式解做比較，這版的效果好很多。

其他的 log10, log2 ，偷用 ln 寫好的公式 (這也是一般人常幹的事)

inline double xlog2(const double x) // slow, 390 / 1266
{ 
    return 1.4426950408889634*xln(x);
}
inline double xlog10(const double x)  // slow, 359 / 1297
{
    return 4.3429448190325176e-1*xln(x);
}