標題猶豫了很久,本想以 [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 問題,整個系列文有興趣可看連結。
綜合以上幾點,採用變數替換之公式,程式碼如下所述。
- typedef unsigned long long u64;
- const double X_LN2 = 0.693147180559945;
- const u64 X_NAN_HEX = 0xfff0000000000001ULL; // not a number
- double xlog(double x, double TOL)
- {
- if(x>0) {
- const u64 hex = *(u64*)&x;
- const u64 t2 = (hex & 0x000fffffffffffffULL) | 0x3ff0000000000000;
- // x = mantissa * 2^exponment
- const double mantissa = *(double*)&t2;
- const int exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
- const double y = (mantissa-1.0) / (mantissa+1.0);
- const double y2 = y*y;
- double yn = y2, cnt=3.0;
- double ans1, ans2=1.0;
- // find nature log (mantissa)
- do{
- ans1 = ans2;
- ans2 += yn / cnt;
- yn*=y2, cnt+=2.0;
- }while(ans2>TOL+ans1 || ans1>TOL+ans2);
- // log(x) = log(mantissa * 2^exponment)
- // = log(mantissa) + exponment * log(2)
- return exponment * X_LN2 + 2.0*y*ans2;
- } else {
- return *(double*)(&X_NAN_HEX);
- }
- }
但注意個小地方,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 ,效果應該會比較快吧?基於此概念程式碼如下
- double xlog2(double x, double TOL)
- {
- if(x>0) {
- const u64 hex = *(u64*)&x;
- const u64 t2 = (hex & 0x000fffffffffffffULL) | 0x3ff0000000000000;
- // (!) x = mantissa * 2^exponment
- const double mantissa = *(double*)&t2 * 0.5;
- const int exponment = (int)((hex >> 52) & 0x7ffULL) - 1022;
- const double y = (mantissa-1.0) / (mantissa+1.0);
- const double y2 = y*y;
- double yn = y2, cnt=3.0;
- double ans1, ans2=1.0;
- // find nature log (mantissa)
- do{
- ans1 = ans2;
- ans2 += yn / cnt;
- yn*=y2, cnt+=2.0;
- }while(ans2>TOL+ans1 || ans1>TOL+ans2);
- return exponment * X_LN2 + 2.0*y*ans2;
- } else {
- return *(double*)(&X_NAN_HEX);
- }
- }
很遺憾的是,這種作法事實上會比較慢,甚至若乘 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 ,
也就是實際上我們可以直接展開前三項,不用再做其它判斷。
- double xlog_3iterm_2(double x)
- {
- if(x>0) {
- u64 hex = *(u64*)&x;
- u64 t2 = (hex & 0x000fffffffffffffULL) | 0x3ff0000000000000;
- // (!) x = mantissa * 2^exponment
- const double mantissa = *(double*)&t2 * 0.5;
- const int exponment = (int)((hex >> 52) & 0x7ffULL) - 1022;
- const double y = (mantissa-1.0) / (mantissa+1.0);
- const double y2 = y*y;
- const double y4 = y2*y2;
- return 2.0*y*(1.0+0.3333333333333333*y2+0.2*y4)+ exponment * X_LN2;
- } else {
- return *(double*)(&X_NAN_HEX);
- }
- }
上面這種展開顯得較差些,這裡使用 Horner Rule 有兩個好處:乘法次數少,重點是浮點數的誤差也少。
- double xlog_3iterm(double x)
- {
- if(x>0) {
- const u64 hex = *(u64*)&x;
- const u64 t2 = (hex & 0x000fffffffffffffULL) | 0x3ff0000000000000;
- // (!) x = mantissa * 2^exponment
- const double mantissa = *(double*)&t2;
- const int exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
- const double y = (mantissa-1.0) / (mantissa+1.0);
- const double y2 = y*y;
- return 2.0*y*(1.0+0.3333333333333333*y2*(1.0+0.6*y2))
- + exponment * X_LN2;
- } else {
- return *(double*)(&X_NAN_HEX);
- }
- }
若要達到 float 之精度 (2^-23 ~= 1.19e-7) ,用一樣的方式去計算,最多需計算 6 次,誤差約為 1.07e-7,一樣套用 Horner Rule 。
- double xlog_6iterm(double x)
- {
- if(x>0) {
- const u64 hex = *(u64*)&x;
- const u64 t2 = (hex & 0x000fffffffffffffULL) | 0x3ff0000000000000;
- // (!) x = mantissa * 2^exponment
- const double mantissa = *(double*)&t2;
- const int exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
- const double y = (mantissa-1.0) / (mantissa+1.0);
- const double y2 = y*y;
- return exponment * X_LN2 + 2.0*y*(
- 1.0 + 0.3333333333333333*y2 * ( // 1/3
- 1.0 + 0.6 * y2 * ( // 3/5
- 1.0 + 0.7142857142857143*y2 * ( // 5/7
- 1.0 + 0.7777777777777778*y2 * ( // 7/9
- 1.0 + 0.8181818181818182*y2 ))))); // 9/11
- } else {
- return *(double*)(&X_NAN_HEX);
- }
- }
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 程式碼約如下。
- double agm(double x, double y, double TOL)
- {
- double an = 0.5 * (x+y);
- double gn = sqrt(x*y);
- double a1, g1;
- do{
- a1 = an, g1 = gn;
- an = 0.5 * (a1+g1);
- gn = sqrt (a1*g1);
- }while(fabs(an-gn)>TOL);
- return gn;
- }
現在所有符號公式都知道了,但,該怎做?
其實沒用到太多技巧,直接對 x 做高精度近似公式,為方便起見,假設 TOL=1e-9,
此時 s = x*2m 裡之 m 至少需為 20 方可滿足此需求, 但別忘了,
這動作實際上是對 x 之浮點數表示法中,exopnment 欄位 + m 而已,
實際上沒必要真的做 * 2m。
- double xlog_high2(double x)
- {
- if(x>0) {
- // 欄位分析
- const double m_mul_ln2 = 13.8629436111989062; // 20 * ln2
- const double PI_D2 = 1.5707963267948966;// PI/2
- // 做高精度近似計算
- const u64 s_hex = *(u64*)&x + 0x0140000000000000; // exp+20, m=20
- const double s = *(double*)&s_hex;
- return PI_D2/agm(1.0, 4.0/s,1e-12)-m_mul_ln2;
- } else {
- return *(double*)(&X_NAN_HEX);
- }
- }
注意,若事先對 x 做 ieee754 欄位切割,再做高精度近似,效果並沒比較好。
- double xlog_high(double x)
- {
- if(x>0) {
- // 欄位分析
- const double PI_D2 = 1.5707963267948966;// PI/2
- const u64 hex = *(u64*)&x;
- const u64 t2 = (hex & 0x000fffffffffffffULL) | 0x3ff0000000000000;
- const double mantissa = *(double*)&t2;
- const int exponment = (int)((hex >> 52) & 0x7ffULL) - 1023;
- // 做高精度近似計算
- const u64 s_hex = (hex & 0x000fffffffffffffULL) | 0x4130000000000000;
- const double m = 20.0;
- const double s = *(double*)&s_hex;
- // const double m_mul_ln2 = 13.8629436111989062; // 20 * ln2
- // return exponment * X_LN2 + PI_D2 / agm(1.0, 4.0/s) - m*X_LN2;
- return PI_D2/agm(1.0, 4.0/s,1e-12) + X_LN2 * ( exponment - m);
- } else {
- return *(double*)(&X_NAN_HEX);
- }
- }
vc 2010, release mode, 一份測時數據供參考 < 別太重視快了多少之類的 >
- log : td = 1609 < math.h >
- xlog : td = 2406 <欄位分析 , 變數替代>
- xlog2 : td = 2422 <欄位分析 , 做移位處理>
- xlog_3iterm : td = 1672 <欄位分析 , 前三項展開 1.4e-4, 使用 Horner Rule>
- xlog_3iterm_2 : td = 1703 <欄位分析 , 前三項展開, 不使用 Horner Rule>
- xlog_6iterm : td = 2078 <欄位分析 , 前六項展開 1.2e-7, 使用 Horner Rule>
- xlog_high : td = 6125 <欄位分析 , 高精度近似>
- 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) 會小於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 欄位拆解做為一些函式技巧,關於其他技巧可能較為深入,甚至可能最後必須找出魔術數字以期加速,這些顯得較為艱澀 ( 找到魔術數字的通常都已發表論文 ),敘述至此結束。
留言列表