標題猶豫了很久,本想以 [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*2裡之 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 * 2, 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) 會小於1,方便起見,用倒數關係 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 欄位拆解做為一些函式技巧,關於其他技巧可能較為深入,甚至可能最後必須找出魔術數字以期加速,這些顯得較為艱澀 ( 找到魔術數字的通常都已發表論文 ),敘述至此結束。

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