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

edisonx

[註 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 = 1010(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 很多都已用進去，然後再進行優化過。

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