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);
}

 

實際上效果和 math.h 顯差很多,原因是這裡的寫法很單純、直接方式去做,

目前查到較有水準的 math.h source code,在不調用 inline asm、用 C 寫之情況下,

大致都是直接在 floating 欄位上進行操作。

 

若對 math.h 有興趣,可上網 google c_prog2.rar ,可得到一些不錯的結果,

這並非唯一的 source ,但可給想設計 math library 之網友一些參考。

arrow
arrow
    全站熱搜

    edisonx 發表在 痞客邦 留言(0) 人氣()