close

標題猶豫了很久，本想以 [C語言數值分析] cmath / math.h 實作提要紀錄 (II) 為標題，

不過這篇真正的重點，筆者本意是較偏向以 ieee754 方式分析 library function，

只是以 log (nature log, ln) 為例做為主要分析對象。

 

所謂的「加速」math library 指的不是加速「內建的」，而是怎麼透過浮點數欄位分析，

去改善自己寫的 math library。

 

在 [C語言數值分析] cmath / math.h 實作提要紀錄 的時候筆者大致將各函式可行性做法推過一次，

這篇做一些補充與加速技巧。

而這裡的「加速技巧」，是可基於 IEEE754 分析方式之函式 ( 所以三角函式第一個被 fired )，

不難，但實用。

 

相容性問題

一些假設先提出

(1) unsigned long long 為 64 bits。

(2) double 符合 IEEE754 64-bit 規範 <剛好手邊的 vc/gcc 大致都符合 > 。

(3) 不考慮特殊數值，如 NAN, INF, NON-NORMALIZE 等問題。

(4) 敘述時，將假設輸入引數 > 0 。

(5) TOL ( 最小容許誤差) 沒特別指明時，是設 1e-12，所以與內建之 math.h 預設 DBL_EPS 有點差異，測時之效能確實也有失客觀性。

另關於筆者一些文章，在「講究標準」之讀者看來，是 garbge，原因是太多假設是 C99 standard 裡沒明文規定 ( 像是指標轉換只能先轉 void* ，再轉到 Type* 之類的)，這部份筆者不多著筆墨。

程式碼以 nature log 為主要示例，其他的僅以文字敘述。

 

難度高之 math.h library

 

個人認為要做到 正確、效率好 的 math library，都與浮點數操作相關(fabs 例外)，諸如

frexp、ldexp、modf、fmod、floor、ceil ，

外加筆者認為大多 compiler 內容會多寫的 is_int，

同時許多 math function 也與這幾支函式之效能有不小之關聯性。

這 7 個浮點數函式大致上不怎麼可能用 C 去做，

不是做不到，用 pointer + bitwise 方式可完成，

差一點用 union + struct 可完成 < 效率會差一點 > ，

但這七個函式目前筆者以 C 方式完成，效能差非常非常多，

推測大致上應都是以組語方式完成，組語應也較為合適，

故這裡不再對其操作做解說。

 

 

Nature Log - 展開式

 

在 wiki 網頁 上， log 所使用的展開式約有四種

 

(1) 直接展開

log(x) = (x-1) - (x-1)2 / 2 + (x-1)3 / 3 - (x-1)4 / 4 + ...

條件為 |x-1| < 1 

(2) Euler Transform

log( x / (x-1) ) = 1 / x + 1 / (2x2) + 1/(3x3) + ....

(3) 變數替代

這公式實際上是推導來的，推導不難，寫起來較麻煩，這裡就不推了。

求 log(x) , 令 y = (x-1) / (x+1), 所以 y 一定會介於 [0,1) 之間。

log(x) = 2 * ( y + y3/3 + y5 / 5 + y7 / 7 +.... )

(4) 高精度近似

log(x) ~=  pi / 2M(1, 4/s) - m log(2)

wiki 上這公式是從牛頓法推導出來的，但其中的 M(1,4/s) 不斷地呼叫根號，實作上顯得較麻煩或不易快。

 

這裡不討論 Euler Transform 與 高精度近似 法，只以「直接展開」與「變數替代」做為說明。有個特性要先說，欲達到相同精度，當輸入引數 x 值愈大時，計算次數就愈多。

 

Nature Log - 基本正確作法


直接展開

log(x) = (x-1) - (x-1)2 / 2 + (x-1)3 / 3 - (x-1)4 / 4 + ...

它有個限制，引數 x 必須滿足 |x-1| < 1 之條件，但大多情況可能都不符合，所以不做任何轉換基本上是辦不到的。

此時可善用 frexp 可達到這部份需求，frexp 主要將引數 x 轉換為

x = mantissa * 2exponment

且用 frexp 切割有個好處是，它返回的小數部份保證在 [0.5, 1) 之間，這裡我們簡寫成以下型式

x = m * 2e

切割好後就是數學問題了

log(x) = log(m*2e) = log(m) + log(2e) = log(m) + e*log(2)

上面的 log(2) 本身是常數，可以藉由事先計算取得，所以實際上真正要算的，只有 log(m) 而已。這是最低限度的作法，不示例。


變數替換

求 log(x) , 令 y = (x-1) / (x+1), 所以 y 一定會介於 [0,1) 之間。

log(x) = 2 * ( y + y3/3 + y5 / 5 + y7 / 7 +.... )

在這裡，雖然替換後可以直接展開，但通常不這麼做，原因有二，一是誤差非常大，另一原因是當 x 較大時，它的收斂速度非常慢，所以事先工作也是和上一個方法一樣，先將 x 分解成 x = m * 2e，接下來的步驟就沒什麼稀奇了。

注意的是，這個公式的收斂速度比上一個「直接展開」的收斂速度快很多！

 

Nature Log - 從 IEEE754 format 分析

 

呼叫 frexp 的目的，一方面確保小數部份滿足公式之條件，一方面可減少迭代次數，

但自己寫的 function 額外呼叫 frexp 之時間成本顯得不少，於是從 ieee754 format 下手。

關於 ieee754 問題，整個系列文有興趣可看連結。

綜合以上幾點，採用變數替換之公式，程式碼如下所述。

 

Code Snippet
  1. typedef unsigned long long u64;
  2. const double X_LN2     = 0.693147180559945;
  3. const u64    X_NAN_HEX = 0xfff0000000000001ULL; // not a number
  4. double xlog(double x, double TOL)
  5. {
  6.     if(x>0) {
  7.         const u64 hex  = *(u64*)&x;
  8.         const u64 t2   = (hex & 0x000fffffffffffffULL)  | 0x3ff0000000000000;
  9.  
  10.         // x = mantissa * 2^exponment
  11.         const double mantissa  = *(double*)&t2;
  12.         const int    exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
  13.  
  14.         const double y = (mantissa-1.0) / (mantissa+1.0);
  15.         const double y2 = y*y;
  16.  
  17.         double yn = y2, cnt=3.0;
  18.         double ans1, ans2=1.0;
  19.  
  20.         // find nature log (mantissa)
  21.         do{
  22.             ans1 = ans2;
  23.             ans2 += yn / cnt;
  24.             yn*=y2, cnt+=2.0;
  25.         }while(ans2>TOL+ans1 || ans1>TOL+ans2);
  26.         
  27.         // log(x) = log(mantissa * 2^exponment)
  28.         // = log(mantissa) + exponment * log(2)
  29.         return exponment * X_LN2 + 2.0*y*ans2;
  30.     } else {
  31.         return *(double*)(&X_NAN_HEX);
  32.     }
  33. }

 

但注意個小地方，ieee754 在存入浮點數時，是以

+/-  1.mantissa * 2exponment

方式存入，用 frexp 時取出之小數保證在 [0.5, 1) 以內，但直接從 ieee754 欄位分析取出時，

它只能保證取出之小數在 [1.0, 2.0) ，故迭代次數肯定比單純用 frexp 還多次；

二擇一的地方在於，呼叫 frexp 的時間，並不比 自行切欄位 + 迭代成本 來得划算。

相對的，如果自己能寫出執行速度快的 xfrexp, 使得取出小數保證在 (0.0, 0.5]，效率會更高，

這部份筆者還沒有心得，本質上都是純粹切欄位進行。

 

 Nature Log - 別進行欄位調整


一種想法是，既然取出之浮點數 mantissa 介於 [1.0, 2.0) 之間，那直接將這個浮點數除 2，

exponment+1 ，效果應該會比較快吧？基於此概念程式碼如下


Code Snippet
  1. double xlog2(double x, double TOL)
  2. {
  3.     if(x>0) {
  4.         const u64 hex  = *(u64*)&x;
  5.         const u64 t2   = (hex & 0x000fffffffffffffULL)  | 0x3ff0000000000000;
  6.         // (!) x = mantissa * 2^exponment
  7.         const double mantissa  = *(double*)&t2 * 0.5;
  8.         const int    exponment = (int)((hex >> 52) & 0x7ffULL) - 1022;
  9.  
  10.         const double y = (mantissa-1.0) / (mantissa+1.0);
  11.         const double y2 = y*y;
  12.  
  13.         double yn = y2, cnt=3.0;
  14.         double ans1, ans2=1.0;
  15.  
  16.         // find nature log (mantissa)
  17.         do{
  18.             ans1 = ans2;
  19.             ans2 += yn / cnt;
  20.             yn*=y2, cnt+=2.0;
  21.         }while(ans2>TOL+ans1 || ans1>TOL+ans2);
  22.         return exponment * X_LN2 + 2.0*y*ans2;
  23.     } else {
  24.         return *(double*)(&X_NAN_HEX);
  25.     }
  26. }


很遺憾的是，這種作法事實上會比較慢，甚至若乘 0.25、exponment+2，效果其實會慢得更明顯。另一種用 bitwise 完成此作法為

(1) 先取出 double x 欄位之 mantissa_x , exponment_x

(2) double m 欄位設定 : mantissa = mantissa_x , exponment = -1  

(3) int exp 取出 : exp = exponment+1

上述之 +/- 1 , 實際代表的是除 2^n ，原本 mantissa 取出範圍為 [1.0, 2.0) , 除二後變成 [0.5, 1.0) ，但即使經過調整後，算出之結果效果還是比之前的還慢，甚至誤差也來得較大。

 

Nature Log - 較低精度時...


在一些對 math library 精度沒那麼要求之情況下，如只需要達到 1e-3 之情形，

我們是可以直接將 TOL 以 1e-3 代入上述之 xlog，但通常會再進一步分析。

 

由於從浮點數欄位分析時，迭代次數最多的情況，

是發生在 mantissa = 1.9999999999999998 ( 欄位全 1 ) 之情況下 ，

配合 TOL = 1e-12，拿這個數字去跑 xlog ，會發現 xlog 需要 12 次計算。

 

相對的，我們再拿 TOL = 1e-3 與 mantissa = 1.9999999999999998 去跑，

發現了 xlog 只需要 3 次計算便可，而且真正的誤差量約是 1.4 e-4 ，

也就是實際上我們可以直接展開前三項，不用再做其它判斷。

 

Code Snippet
  1. double xlog_3iterm_2(double x)
  2. {
  3.     if(x>0) {
  4.         u64 hex  = *(u64*)&x;
  5.         u64 t2   = (hex & 0x000fffffffffffffULL)  | 0x3ff0000000000000;
  6.         // (!) x = mantissa * 2^exponment
  7.         const double mantissa  = *(double*)&t2 * 0.5;
  8.         const int    exponment = (int)((hex >> 52) & 0x7ffULL) - 1022;
  9.         const double y = (mantissa-1.0) / (mantissa+1.0);
  10.         const double y2 = y*y;
  11.         const double y4 = y2*y2;        
  12.         return 2.0*y*(1.0+0.3333333333333333*y2+0.2*y4)+ exponment * X_LN2;
  13.     } else {
  14.         return *(double*)(&X_NAN_HEX);
  15.     }
  16. }

 

上面這種展開顯得較差些，這裡使用 Horner Rule 有兩個好處：乘法次數少，重點是浮點數的誤差也少。

 

Code Snippet
  1. double xlog_3iterm(double x)
  2. {
  3.     if(x>0) {
  4.         const u64 hex  = *(u64*)&x;
  5.         const u64 t2   = (hex & 0x000fffffffffffffULL)  | 0x3ff0000000000000;
  6.  
  7.         // (!) x = mantissa * 2^exponment
  8.         const double mantissa  = *(double*)&t2;
  9.         const int    exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
  10.  
  11.         const double y = (mantissa-1.0) / (mantissa+1.0);
  12.         const double y2 = y*y;
  13.  
  14.         return 2.0*y*(1.0+0.3333333333333333*y2*(1.0+0.6*y2))
  15.             + exponment * X_LN2;
  16.     } else {
  17.         return *(double*)(&X_NAN_HEX);
  18.     }
  19. }

 

若要達到 float 之精度 (2^-23 ~= 1.19e-7) ，用一樣的方式去計算，最多需計算 6 次，誤差約為 1.07e-7，一樣套用 Horner Rule 。

Code Snippet
  1. double xlog_6iterm(double x)
  2. {
  3.     if(x>0) {
  4.         const u64 hex  = *(u64*)&x;
  5.         const u64 t2   = (hex & 0x000fffffffffffffULL)  | 0x3ff0000000000000;
  6.  
  7.         // (!) x = mantissa * 2^exponment
  8.         const double mantissa  = *(double*)&t2;
  9.         const int    exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
  10.  
  11.         const double y = (mantissa-1.0) / (mantissa+1.0);
  12.         const double y2 = y*y;
  13.  
  14.         return exponment * X_LN2 + 2.0*y*(
  15.             1.0 + 0.3333333333333333*y2 * ( // 1/3
  16.             1.0 + 0.6 * y2 * ( // 3/5
  17.             1.0 + 0.7142857142857143*y2 * ( // 5/7
  18.             1.0 + 0.7777777777777778*y2 * ( // 7/9
  19.             1.0 + 0.8181818181818182*y2 ))))); // 9/11
  20.     } else {
  21.         return *(double*)(&X_NAN_HEX);
  22.     }
  23. }

 

Nature Log - 高精度近似法


公式為 log(x) ~=  pi / 2M(1, 4/s) - m log(2)，下面這些純粹在翻譯 wiki 東西。

有幾個符號要講，s、m、M() function。

先挑  m ，使得 s = x * 2m    ，同時假設有效位數是 p 位數，

 s 必須滿足， s = x*2m > 2p/2 。，一般而言會挑 m=8。

 

若 x=1, m=8 , 代表精度為小數後 p=2*8=16 位，也就是1.5e-5，

若要達 1e-12, p 必須 = 40 ，此時 m 必須挑 20。

 

再舉個例，

若 x = 20 時，m 挑 16，則可得到 s = x*2m = 20*216，

滿足  s = x*2m > 2p/2  之最大 p 為 28，

此時的精度是 2^-28 ~= 3.7e-9。

 

重點是 M() ，中文叫 算術-幾何 平均數，AGM，不停在調用根號，

AGM 程式碼約如下。

 

Code Snippet
  1. double agm(double x, double y, double TOL)
  2. {
  3.     double an = 0.5 * (x+y);
  4.     double gn = sqrt(x*y);
  5.     double a1, g1;
  6.     do{
  7.         a1 = an, g1 = gn;
  8.         an = 0.5 * (a1+g1);
  9.         gn = sqrt (a1*g1);
  10.     }while(fabs(an-gn)>TOL);
  11.     return gn;
  12. }

  

 現在所有符號公式都知道了，但，該怎做？

其實沒用到太多技巧，直接對 x 做高精度近似公式，為方便起見，假設 TOL=1e-9，

此時 s = x*2m 裡之 m 至少需為 20 方可滿足此需求， 但別忘了，

這動作實際上是對 x 之浮點數表示法中，exopnment 欄位 + m 而已，

實際上沒必要真的做 * 2m。

 

Code Snippet
  1. double xlog_high2(double x)
  2. {    
  3.     if(x>0) {
  4.         // 欄位分析
  5.         const double m_mul_ln2 = 13.8629436111989062; // 20 * ln2
  6.         const double PI_D2 = 1.5707963267948966;// PI/2
  7.  
  8.         // 做高精度近似計算
  9.         const u64 s_hex = *(u64*)&x + 0x0140000000000000; // exp+20, m=20
  10.         const double s =  *(double*)&s_hex;
  11.         return PI_D2/agm(1.0, 4.0/s,1e-12)-m_mul_ln2;
  12.     } else {
  13.         return *(double*)(&X_NAN_HEX);
  14.     }
  15. }

 

 注意，若事先對 x 做 ieee754 欄位切割，再做高精度近似，效果並沒比較好。

 

Code Snippet
  1. double xlog_high(double x)
  2. {    
  3.     if(x>0) {
  4.         // 欄位分析
  5.         const double PI_D2 = 1.5707963267948966;// PI/2
  6.         const u64 hex  = *(u64*)&x;
  7.         const u64 t2   = (hex & 0x000fffffffffffffULL)  | 0x3ff0000000000000;
  8.         const double mantissa  = *(double*)&t2;
  9.         const int    exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
  10.  
  11.         // 做高精度近似計算
  12.         const u64    s_hex = (hex & 0x000fffffffffffffULL)  | 0x4130000000000000;
  13.         const double m = 20.0;
  14.         const double s =  *(double*)&s_hex;
  15.         // const double m_mul_ln2 = 13.8629436111989062; // 20 * ln2
  16.         // return exponment * X_LN2 + PI_D2 / agm(1.0, 4.0/s) - m*X_LN2;
  17.         return PI_D2/agm(1.0, 4.0/s,1e-12) + X_LN2 * ( exponment - m);
  18.     } else {
  19.         return *(double*)(&X_NAN_HEX);
  20.     }
  21. }

 

vc 2010, release mode, 一份測時數據供參考 < 別太重視快了多少之類的 >

Code Snippet
  1. log           : td = 1609 < math.h >
  2. xlog          : td = 2406 <欄位分析  , 變數替代>
  3. xlog2         : td = 2422 <欄位分析  , 做移位處理>
  4. xlog_3iterm   : td = 1672 <欄位分析  , 前三項展開 1.4e-4, 使用 Horner Rule>
  5. xlog_3iterm_2 : td = 1703 <欄位分析  , 前三項展開, 不使用 Horner Rule>
  6. xlog_6iterm   : td = 2078 <欄位分析  , 前六項展開 1.2e-7, 使用 Horner Rule>
  7. xlog_high     : td = 6125 <欄位分析  , 高精度近似>
  8. xlog_high2    : td = 6031 <無欄位分析, 高精度近似>

 數據列出來的主要目的是想講，這種分析方式還是沒比 math.h 裡的來得快，比較值得講的是，高精度近似法事實上本身就是「近似」，在精度要求較小之情況較適合用，在這種情形下，agm 裡頭，sqrt 之算法也可做適當之加速，這裡便不再贅述。

一樣，log 會算之後，log10 用換底公式可完成。

 

其他 math.h 算術類技巧

 

這類函式含 exp、pow、sqrt，分析步驟方式都如 log ，

(1) 找適當的展開函式

(2) 可否將 func( mantissa * 2 exponment) 化簡？

(3) 若可化簡，取得 ieee754 欄位，開始進行分析。

 以下只做部份提示。

 

sqrt(x)

 

已知 S , 欲求 x = sqrt(S), 牛頓法推導如下

x2 = S , 令 f(x) = x2-S , 則 f'(x) = 2x

可設起始點為 0.5 * S ，a0 = 0.5*S

an+1 = an + f(an) / f'(an) = an - (an2-S) / (2an)

     = an - 0.5 * (an - S/an)

     = 0.5 * (an + S/an-1)

直到 |an+1 - an| 小於 TOL 停止。

----

sqrt(x),  x 愈大收斂愈慢，故 x 愈小愈好。

令 x = m * 2e ， 1<=m < 2

sqrt(x) = x^0.5 = (m*2e) ^ 0.5 

sqrt(x) = (m)^0.5 * (2e)^0.5 = m0.5 * 2e/2


所以就變成了求 sqrt(m)。

注意事項：必須做 x < 0 之額外判斷。

 

exp(x)

 

taylor series, exp(x) = 1+x+x2+x3+...+xn　，

第 n 項與第 n+1 項誤差小於 TOL 視為收斂，x 愈大收斂愈慢，

甚至會算到 INF。

 

令 x = m * 2e，

exp(x) = exp(m*2e) = exp(m) * exp(2e)

           = exp(m) * e * exp(2)

其中 exp(2) 為常數，分解 m, e 後，只需求 exp(m) 即可。 

注意事項 : 當 x < 0 時，exp(x) 會小於１，方便起見，用倒數關係 exp(x) = 1.0 / exp(-x) 處理掉 x < 0 之情況。

 

Power

 

這幾個裡面，power 是個最麻煩的函式。

xy = exp ( ln (xy) ) = exp (y * ln(x) )

注意到會先以 ln 求 ln(x)，所以 x 本身並不能為 0、不能為負 ，

但考慮到正確性情況時，如 -2.0 4 答案是 16，所以，

必須再寫一個函式去判斷指數部份 y 是否為整數，如果是的話就呼叫 fast power 求解；

如果不是整數、x 又小於 0 ，那這就不會有正確 ( 驗證過，是此規則無誤 )，

下面說的情況是 y 非整數，x 大於 0 之情況

 

xy 我們可以很偷懶的用 exp(y * ln(x) ) 去求，剛好 ln 和 exp 剛剛都講過技巧，所以可以算出來。

 

另一種方式，再將 x, y 從 IEEE754 欄位拆出來，看能不能化簡。

令 x = m1 * 2e1 ，y = m2 * 2e2 ，則

xy = exp( y * ln(x) ) = exp( m2 * 2e2  * ln(m1 * 2e1) )

     = exp(  e1 * m2 * 2e2  * ln(m1 * 2) )

     = exp(  e1 * m2 * 2e2  * ( ln(2) + ln(m1) ) )

令 ln(2) + ln(m1) = z, 

xy = exp(  e1 * m2 * 2e2  * z ) 

     = e2 * exp(e1) * exp(m2) * exp(2) * exp(z)

上面只有 exp(2) 是常數，但如此下來時，等於是間接呼叫了 3 次 exp 函式、1 次 ln 函式，

這與剛剛的 1 次 exp 函式、1次 ln 函式速度差異簡直過大，原因在於，過度化簡。

事實上上述步驟化到 (令 z) 時就不用再化了，也就是

xy = exp(  e1 * m2 * 2e2  * ( ln(2) + ln(m1) ) )

     = exp( e1 * y * ( ln(2) + ln(m1) )。

 

 注意事項 1 : 個人認為大多編譯器會實作 is_int(x) 這裡是其中一原因，當然調用 ceil(x)==x 或 floor(x) == x 也可完成，唯應是 ceil(x)、floor(x) 內部調用 is_int(x) 之可能性較高。

注意事項 2  : 上述推導可知，實際上 pow 函式確實為相當耗時，一般在指數為整數時大多自己寫 fastpower ( 需考慮指數為負數)，也別企圖以 pow(x,0.5) 企圖取代 sqrt(x)

注意事項 3 : 目前筆者不知道有沒有對於 pow function 具更有效率之解法，上述概念實作後約是內建之 2~3 倍慢。


小結

 

本文以 IEEE754 欄位拆解做為一些函式技巧，關於其他技巧可能較為深入，甚至可能最後必須找出魔術數字以期加速，這些顯得較為艱澀 ( 找到魔術數字的通常都已發表論文 )，敘述至此結束。

arrow
arrow
    全站熱搜

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

    E-mail轉寄