微剋多資訊

 找回密碼
 註冊

Login

Login

搜索
回覆 0則 瀏覽 1729篇

[愛迪生系列文] 02 - 談談冪次(pow)

該用戶從未簽到

發表於 2012-12-14 03:54 | 顯示全部樓層 |閱讀模式
本文章最後由 edisonx 於 2012-12-14 04:20 編輯

[註 1] 系列文不會講基本語法與概念,這部份筆者做太多了,有點膩。
[註 2] 系列文都屬原創基於知識無價原則,您可以自由使用、修改,唯基於網路品德,請註明作者與出處,避免誤觸智財權

[註 3] 系列文參考書籍可能很多,唯筆者一些經驗與整理後之發表。


[1] C/C++ 調用之 pow 函式



首先在 C++ 之數學函式庫是 cmath,C 語言之數學函式庫是 math.h,對於 pow 函式,C/C++ 兩者差別在於,C++ 裡面多載了三份 pow

long double pow(long double x,long double y);
float pow(float x, float y);
double pow(double x, double y);

但 C 語言裡面沒有多載的概念,所以只有一份

double pow(double x, double y);

在呼叫調用時,確保資料型態之正確性,與資料型態隱含轉換所帶來之影響,這裡就不再贅述。


[2] 其實.. pow 速度很慢!



常看到一些人要算某個數的平方時

double x =1.2;
double y= pow(x, 2.0);

或是要算根號時

double x =1.2;
double y= pow(x, 0.5);

這速度其實都很慢!原因在本文後面會說明。一般而言,如果要算平方的話其實寫 y = x*x 會快很多!根號就直接寫 y = sqrt(x) ; 也一樣會快很多!注意,是「快很多」,不是「快一點點」。


[3] 幕次為整數時候,還是乾脆點自己寫一個吧!



當 pow(x, y) 裡面的 y 是整數的時候,呼叫自己寫的一定會比較快!無論自己寫的函式有多醜,只要針對幕次是整數寫一個函式出來,就一定比原本的 pow 還要快。注意的是,一般錯誤的寫法是

double powi(double x,int n){
     
double rst=1.0;
     
while(n--)rst*=x;
     
return rst;
}


錯誤的地方在於,n 可能本身為負數,帶進去時會出現錯誤。這裡更正方式其實很簡單,假設 n 為正數, a^(-n) = 1.0 / (a^n),故寫下以下程式碼。

double powi(double x,int n){
     
double rst=1.0;
     
int m = (n>=0) ? n: -n;
     
while(m--)rst*=x;
     
return (n>=0)? rst : 1.0/rst;
}


修正成這樣就正確了。這段程式碼,在幕次為整數的情況,一定比 pow 來得快!有興趣可自己做測時。


[4] 傳說中的 fast_power



傳說中的 fast_power ,是在幕次為整數的時候使用,速度比上面說的 powi 還快!我們以 2^10為例,它的概念如下。

(0)  初始化,exp = 10(dec) = 1010(binary),rst = 1,tmp = base = 2;
(1)  第一次,exp = 101
0(binary),rst 不變 , rst = 1,
       tmp = tmp * tmp = 4, exp>>=1, exp=0101
(2)  第二次,exp = 0101(binary),rst = rst * tmp = 1 * 4 = 4 ,
       tmp = tmp * tmp = 16, exp>>=1, exp=0010
(3)  第三次,exp = 0010(binary),rst 不變 , rst = 4
       tmp = tmp * tmp = 256, exp>>=1, exp=0001
(4)  第四次,exp = 0001(binary),rst = rst * tmp = 4 * 256 = 1024 ,
      tmp = tmp * tmp = 65536, exp>>=1, exp=0000
(5) 第五次,exp 變 0000 ,即 0,結束,最後 rst = 1024 即為答案。

轉成 C 語言如下。

double fast_pow(double x,int n){
     
double rst=1.0;
     
double tmp = x;
     
unsigned m = (n>=0) ? n: -n;
     
while(m){
         
if(m&1)rst*=tmp;
         tmp*=tmp;
         m>>=1;
     }
     
return (n>=0)? rst : 1.0/rst;
}


[5] 我想實作 C 語言的 pow 函式



我們想一下 pow(a, b) ,在 a, b 都是浮點數的時候,這動作是怎麼完成的。首先,pow(a, b) 事實上沒辦法直接算出來,它是間接轉換的。

令 a^b = C ,
兩邊取 ln 得到 --->  ln(a^b) = C ---> b * ln(a) = ln(C),
再對  b * ln(a) = ln(C) 兩邊取 exp,得到
exp( b * ln(a) ) = exp( ln(C) ) ---> exp( b * ln(a) ) = C = a^b


說穿了,當我們在算 a^b 時,實際上是在算 exp(b * ln(a)),
那 exp、ln 函式又怎麼計算的呢?



[6] exp 函式



我們可以在 wiki上找到 exp 之 taylor series,

exp(x) = 1 + x + x^2/2! + x^3 / 3! + ....

然後關鍵不是要計算幾項,而是要計算,當第 n 項與第 n+1 項之誤差小於我們指定的誤差(像是 1E-9,這種誤差稱最小可接受誤差,eps)時,就停止計算。先照著它的定義跑,我們可以寫出下列程式碼

double xExp(const double x,const double eps)
{
     
double ans1 ;
     
double ans2 = 1.0;
     
double fact=1, xn=x, cnt=2.0;
     
do{
         ans1=ans2, ans2 = ans1+xn/fact;
         fact*=cnt;
         xn = xn*
x;
         cnt+=1.0;
     }
while((ans1 > ans2 + eps) || (ans2 > ans1 + eps));
     
return ans2;
}

重點又來了,上面的程式碼其實效能很差,原因在於當 x 愈大的時候,所需要的計算時間就愈長,而且可以算的範圍,比內建的 exp 還要小很多!問題出在哪?

問題出在,其實 math.h  / cmath 之 exp ,有用公式轉換過,我們將浮點數 floating 拆成整數部份 p1 與小數部份 p2,


floating = p1 + p2,
exp(floating) = exp(p1 + p2)
exp(floating) = exp(p1) * exp(p2)


然後 exp(p1) 部份,因為 p1 是整數,所以可以用 fast_power 快速求出其值,而 p2 是小數,再代入上面的 xExp ,速度會較快。故程式碼就變成下面。

double xExpFast(const double exponment,const double eps)
{
     const double Euler = 2.718281828459045;
     double rst=1.0;
     double p1=ceil(exponment), p2=exponment-p1;
     if(exponment> 709.0){
         p1=1.0, p2=0.0;
         returnp1/p2;
     } else if(exponment< -709.0) {
         return0.0;
     } else {
         return xExp(p2, eps) * fast_pow(Euler, (int)p1);
     }   
}


[7] Natrue Log 函式



exp 求完了,那 nature log 呢?ln (x) 我們可以這麼算

(1) 令 t = (x-1) / (x+1)
(2) 則 ln(x) = 2t * (  1 + t^2/3  +  t^4 / 5  + t^6  / 7 + .... )


一樣,照著上面的定義有很多數值求出來會有問題,因實際上 library 裡之 nature log 有用過公式轉換過一層,不然很多數值會求不出來。假設要求的是 ln(x),且 x = A * 2^B,接著

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

ln(2) 是常數可以事先寫好,而這時候 ln(A) 再套入上述之公式即可。那關鍵是,怎麼將 x = A * 2^B ?? 事實上 math.h 裡面準備了一個函式等著了,叫 frexp。

frexp 這個函式是從浮點數欄位(很有可能是  ieee754) 裡面分析拿到的,而且這個函式很不建議自己做,原因是自己做速度超慢,另一個原因是做這類型函式等於是自己再定義系統的浮點數表示法。無論如下, log 函式如下所示。


double xln(const double x,const double eps){   
     
int B;
     
constdouble A = frexp(x, &B); // x = A* 2^B

     
// findln(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*6.9314718055994529E-1;
}



[8] 回歸到 pow



如上面所提到的,pow(a, b) 事實上是呼叫 exp(b * ln(a)),所以再套入即可。

double xpow(const double base,const double exponment,constdouble eps)
{
     return xExp( exponment* xln(base, eps),eps);
}

上面這部程式其實還有一點問題,想一下 (-2)^3 結果是多少?用 pow(-2,3) 其實不同 compiler 答案也是不一樣,一種答案是 -8 ,一種答案是 INF 之類的東西,原因在於 (-2)^3 最後會被化成 exp( 3 * ln(-2)),關鍵在 ln(-2) 的部份會是無定義行為,所以如果要真的正確的話,必須還要再額外判斷,當 base 是負數、而 exponment 為整數時,這時候就該調用 fast_power 答案才會正確,不能直接用 exp( b * ln(a) ) 完成。這裡就不再示範



[9] 補充說明



(1) 這篇文章寫的 pow 是一種實作方式,使用的是 taylor series ,另一種是以 Pade approache 方式展開,這裡比較深入,就不再提了。

(2) 除了指數是整數的情況下,不要期待會比內建的 pow  還快。

(3) 其他的數學函式沒意外的話,不論是精度還是速度,都比自己寫的還快!包含在網路上找到的什麼 fast_sin、fast_sqrt 、fast_sqrt_inv 等,這些技巧 library 很多都已用進去,然後再進行優化過。


愛迪生系列文第二篇就至此,有興趣討論的版友歡迎留言討論。



您需要登入後才可以回帖 登入 | 註冊

本版積分規則

小黑屋|Archiver|微剋多資訊(MicroDuo)

GMT+8, 2016-12-6 18:12

Discuz! X

© 2009-2016 Microduo

快速回覆 返回頂部 返回列表