[活動]《我的包包,裝我的旅行》今年暑假,送你去度假![公告] 豐掌櫃《超級吸金王》活動賽況排行及得獎名單 (7/19第二批名單公佈)[公告] 即日起,MIB也能變購物金囉![公告] 豐掌櫃《最佳銷售王》賽況排行 (6/15 得獎名單公佈)[公告] 痞客邦 PIXNET MIB(MONEY IN BLOG)部落格廣告分潤計劃申請流程調整

二年前瘋狂研究 math.h 裡之各 function 如何實作,

一般見得了人之 function 為求速度,都直接與 IEEE754 format 做為操作基準,

這裡提的是「見不了人」的,也就是速度上有所懸殊之方式 < 但肯定比學校作業要求好多了 ...>

 

猜我大概也沒什麼時間精力完成這系列文章之所有實作,留個筆記底,

避免日後 < 牛奶事件 >[文末註] 發生時還能有所參考。

 

切記一點:Euler’s Infinite Series ,此系列之 series 用的都是連乘積,

在理論上收斂速度會比 Taylor series 還要快,但不代表它的速度會比 Taylor series 還快,

原因是它在計算之過程中,所費的成本比 Taylor series 還多。

 

 

壹、浮點數運算類

 

1. ceil 

 

(1) 抓 ieee754 - exp,轉回實際值。

(2) exp < 0,直接傳回 1 - signed bit ( or NOT signed bit)

(3) exp > 0,先將 mantissa 左移 exp 個 bits,判斷左移後是否為零。為零:直接傳回原始值;不為零:移出部份+1,沒移出部份清0。

 

2. floor 

 

分析方式與 ceil 同。

 

3. fabs 

 

可用 pointer 或 if-else,建議以 if-else 方式實作即可,速度未必較慢。

 

4. fmod 、modf  

 

 直接去看算盤本裡面怎麼處理 IEEE754 浮點數除法問題比較實際。

 

5. frexp、ldexp

 

這二個函式只是在考驗對 IEEE754 熟悉度,與指標、bitwise 之熟練度。

 

 

貮、三角函式類


 

 

所有的 fast_sin / fast_cos 誤差都很大,且不會比 cmath / math.h 之 sin/cos 來得快。

 

1. sin

 

taylor series 展開, sin (x ) = x - x^3/3! + x^5/5! - .....,

先將 x 轉到主幅角後 (mod 2PI) 再做運算收斂較快,直到迭代誤差小於 eps 為止,

建表法可以做,1024 份,最大誤差約在 0.1 還是 0.01 < 忘了 >。

 

另 sin 雖 2 倍角公式不能用 ( sin2x = 2sinxcosx ),

但幸運的是 sin3x , sin5x 都可以拿來用,x 範圍縮小後,

自然可增加收斂速度,但切記使用倍角數愈大,

所產生之誤差很可能就愈大。

sin(3x) = 3sinx - 4sinx^3  = sinx * (3 - 4 sinx)

令 sin(x) = y

sin(5x) = 16*y^5 - 20*y^3 + 5*y = y*(5+y^2*(-20+16*y^2))

輔助特性: x 接近 0 時, sin(x) = x 。

 

2. cos

 

與 sin 一樣,用 taylor series 展開, cos (x) = 1 - x^2/2! + x^4/4! - ....,

先將 x 轉到主幅角後 (mod 2PI) 再做運算收斂較快,直到迭代誤差小於 eps 為止,

另 cos 可利用二倍角公式: cos 2x = 1 - cosx ^ 2 ,再加快收斂速度 < 精度會少一些 >,

輔助特性:x 接近 0 時,cos(x) = 1 - x。

 

* 另外建議寫一個 function,可以在同一 loop 裡計算 sin ( x ) 與 cos ( x ),

   在圖形學裡面通常這兩個三角函數值會同時用到。

 

3. tan

 

tangent 用 taylor 展開會很麻煩,直接在同一個 loop 裡面算 sin (x) 與 cos(x),

做相除動作,另外還可拿一份半角公式出來用:

tan(2x) = sinx  /  (  1 + cosx) = ( 1 - cosx) / sinx,

算出來精度會降低,但速度應會快一點點

 < 不會快太多,因會增加 brach - prediction 負擔 >

 

4. cot 、sec、csc *

 

這三個 library 裡面沒有,直接做倒數運算即可。

 

 

參、反三角函式類

 

 

反三角函式類大多是基於 arcsin 或 arctan 做為轉換基準,意思是可以發展 arcsin,

也可以發展 arctan,再做其他轉換,筆者習慣以 arctan 為基準。

 

1.  atan

 

使用 taylor series 展開,x - x^3/3 + x^5 / 5 - x^7 / 7 +...

唯需注意一點,這份 taylor 只適用於 |x| < 1 之情況下,

大於 1 時必須使用這特性: tan(x) = PI/2 -  tan(1.0 / x)

愈接近1,它的收斂速度就愈慢,

值得注意的是,在接近於 1 時,反而以 euler series 收斂速度較快。

另外半角公式問題,由於有用到根號,故實際上不會拿來做應用。

輔助特性:0 < x < 0.00000762939453,tan(x) = x。

 

 

2. asin

 

asin 可以用 taylor series 直接展開,

arcsin(x) = x + (x^3/3)*(1/2) + (x^5/5)*(1*3/2*4) + (x^7/7)*(1*3*5/2*4*6) + ...

 

3. acos

 

先算 y = asin(x), acos(x) = PI/2 - y

 

4. atan2

 

atan2(x, y),實際上就是 atan( y / x )。

 

  

 

肆、一般算術類

 


1. exp

 

這問題我覺得學校的範例真的是害死人,只講一半不到的東西。

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

(0) 事先驗證: x > 709 傳回 1.0 / 0.0 , x < -709 傳回 0.0 

(1) 如果 x < 0 的話,exp(x) = 1.0 / exp(-x),拿這點做轉換,直接處理正數部份就好。

(2) 先將 x 切成成整數部份 (A) 與小數部份 (B), x = A + B, 

      exp(x)  = exp(A+B) = exp(A) * exp(B),分二次算會比較快。

(3) exp(A) ,由於 A 是整數,可以做 fast power ( 2.71828.... ^ A 的 fast power)

(4) 再將 exp(B) 套入 taylor series ,收斂到 eps,算出來結果再進行 exp(A) * exp(B) 答案即為所求。

 

注意,上面一定要這樣拆解,如果不拆解的話算不到 exp(70x) 這麼大,同時計算之浮點數誤差也很大。

事實上 library 裡面應也是用這種作法 ( 驗證過誤差,一樣)

 

2. log

 

自然對數,也是用 taylor series 展開,

注意它有兩種表達方式

ln(z) = (z-1) - (z-1)^2/2 + (z-1)^3/3 -..., for 0<z<2 (a)

ln(z) = 2 * sum( 1/(2n+1) * [(z-1)/(z+1)]^(2n+1)] ) , for integer n, z>0 (b)

(a) 應會較慢一點,但 (b) 只比 (a) 好一點,必須要用一點技巧化開。

frexp 函式主要是將一個浮點數 x 化成 x = A * 2^B  (不知道的話請上 c++ reference 查),

所以可以事先做這事

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

其中 ln(2) 已經是常數了,ln(2) =  6.9314718055994529E-1; 

問題就變成解決 ln(A) ,收斂速度快許多。

 

3. log10

 

以10為底的對數,直接基於 log10(x) = ln(x) / ln(10) 做對換。

 

4. log2 * 

以 2 為底的對數,在標準函式庫裡面沒有,也是基於 log2(x) = ln(x) / ln(2) 做對換。

若只是要取整數部份的話,可看這份文章,推用 debruijn。

 

5. sqrt

 

直接使用巴比倫法(可視為牛頓法的特例),

利用 x(n+1) = 0.5 * [ x(n) + S / x(n) ] 迭代求解,小於 eps 視為收斂 ,

根據算術平均 >= 幾何平均,初始值可以給 x / 2。

codeproject 也有人寫了13份 sqrt 之原始碼

 

6. fast_pow *

library 裡從沒有 fast_pow 這東西。

fast_pow 和 pow 最大差別在於,fast_pow 是拿來算,指數為整數次方的,

速度差非常多。大概如下

 

double fast_power(double base, int exponment)
{
    int e = exponment >= 0 ? exponment : -exponment;
    double tmp = base, rst=1.0;
    while(e){
         if(e&1) rst*=tmp;
         tmp*=tmp;
         e>>=1;
    }
    return (exponment >=0) ? rst : 1.0 / rst;
}

 

7. pow

 

pow 函式是我寫過效率最差的內建函式,猜測應有其他方式去寫。

我們以 a^b 為例,首先 b 可拆在整數(X)+浮點數(Y),變成 b = (X+Y),

a^b = a^(X+Y) = a^X * a^Y,

其中 a^X 可以用 fast_pow 去算,重點在 a^Y。

利用對數換底公式,可以得到,

a^Y = exp ( ln ( a^Y ) ) = exp ( Y * ln (a) ),

裡面用到了 exponment function 與 nature log function .

這兩個 function 我們已在上面提過,故可順利解出。

 


伍、雙曲函式

 

 

這裡的分析和三角函式很像,採 taylor series,

也有人拉 sinh(x) = 0.5 * (exp(x) - exp(-x)) 等定義做,但需一點技巧。

 

 1. sinh

 

sinh(x)  = x + x^3/3! + x^5/5! + ... , for all x

誤差到 eps 視為收斂。

 

另一組為 sinh(x) = [ e^2x -1 ] / 2*e^x,

而 e^2x = (e^x)* (e^x),問題轉成了求 e^x,

速度比上述方法慢一點。

 

2. cosh

 

cosh(x) = 1 + x^2 / 2! + x^ 4 / 4! + ...., for all x

誤差到 eps 視為收斂,另外 cosh 應善用其二倍角公式加快收斂

cosh(2x) = 1+2*sinh(x)^2 = 2 * cosh(x)^2 - 1。

但 sinh(2x) = 2*sinh(x) * cosh(x),故 sinh 二倍角助益不大。

 

另一組為 cosh(x) = [ e^2x +1 ] / 2*e^x,

而 e^2x = (e^x)* (e^x),問題轉成了求 e^x。

速度比上述方法慢一點。

 

3. tanh

 

tanh 之 taylor series 使用到 Bermoulli number 較為麻煩,

但也沒必用同時求 sinh(x), cosh(x),直接化成

tanh(x) = [ e^2x - 1 ] / [e^2x + 1] 方式來做較佳。

 

陸、其他函式

 

上述之說明應已將 cmath / math.h 裡之函式跑過一遍,親自實作後會發現有幾個重點

 

1.  查函式定義、定義域、值域。

2.  查有沒有相關之 series (目前一般最常查的的確是 taylor series,但不代表它是最好的,只是最方便) 。

3.  了解該函式之數學轉換 (如商數關係、倍角公式等)。

4.  了解該函式之數值圖形,哪部份的切線斜率較大,就代表那裡的收斂速度快。

 

至於 coth、arcsinh 等函式,由於分析步驟都大同小異,這篇便不再贅述。

 

 

 

 

 

 

[註] 所謂 < 牛奶事件 >,指的是筆者 nb 在論文第四章完成時,被一杯牛奶從 keyboard 那裡整個倒翻,整個 nb 報銷,一切從零開始。

Posted by edisonx at 痞客邦 PIXNET Guestbook(2) 引用(0) 人氣()


open trackbacks list Trackbacks (0)

留言列表 (2)

Post Comment
  • novus
  • fsin fcos 從蠻久以前就已經是 x86 內建指令了。

    查表是有技巧的,不用很大的表就可以達到 20 bit 以上的精度,已經和 float 的精度相近。

    對於不需要高精度的程式,查表通常比 fsin fcos 指令還快很多,這是 x86 有趣的地方--患了功能過多症候群的硬體指令往往比軟體慢....
  • fsin / fcos 讓我想追問幾個問題和 dll 相關的問題 < 有點跳 tone 的話請見諒 > 。 以手邊幾個 compiler 例子來講, vc 裡之 sin 應是放在 MSCVxxx.dll ,gcc 自己又有一份 glibc ,裡面也有 math.h 開放之實作,另外 C:\Windows\System32\ntdll.dll ,dumpbin.exe (VC tool)導出後裡面也有 sin,我能推測的想法是, vc 用的 math.h 自己寫好放在 mscvxxx.dll 裡;gcc 用的是自己寫的 glibc,那 ntdll.dll 裡面的 sin 就是組語用的 fsin 嗎(只是再包成 dll 讓以後會用到的人用而已?)?若它真的是快的,為何各家 compiler 還要自己再弄一份出來?(我能想到的原因是:別人包出來的會比較快,或比較安全)

    查表問題是我最大納悶之處,之前實作一份 [1024+1] 大小之表,
    http://edisonx.pixnet.net/blog/post/83786128 最後一個 section,建起 [0, 2PI]之表,測精度大概在 6E-3,用的是一次線性差值法做修正,且表的大小再上去,能修正的精度非常有限。目前想到可以改善的大概只剩從 [0,2PI] 再轉到 [0, PI/2],但轉過去後又怕在正負號判別時效能被 branch prediction 吃掉,所以猜測您指的可能與此法有所差異,不知這裡是否有哪些 keyword 或文章可供參考?

    最後謝謝 novus 的不吝指導,感謝。

    edisonx replied in 2012/04/17 23:08

  • novus
  • ntdll.dll 裡面的有可能是靜態連結進去的,也有可能是額外實作,真正的原因我不知道,總之都不是太奇怪的事。因為msvcrt「理論上」不屬於windows系統的一部分,win2000 以前的電腦甚至必須額外安裝(通常安裝VS、Offic或其他軟體時順便),ntdll.dll總不能依賴msvcrt才能運作。至於使用靜態連結的理由可能是安全,或者避開DLL hell,我也不曉得。

    Windows上的gcc分兩派,MinGW派的人還是使用msvcrt;而 Cygwin 派的人則堅持要自己弄一套GNU的。至於Linux上則更加多元,C Lib一大堆,也沒有誰好誰壞的問題,有標準作法的東西大家做起來都差不多。

    查表的技巧 Hint:假如你有一套0~360間隔20度的表,和一套0~20度間隔1度的表,總共才多少筆資料?想一下用合角公式你可以回答多少角度sin、cos的問題?

    (如果我每次在網路上回答問題都可以收 $100,那麼就不怕物價上漲了。)
  • 這下真的豁然開朗了!!真是感謝 novus 賜教 !!

    edisonx replied in 2012/04/18 01:01

You haven’t logged in yet, please use guest status to leave message. You can also log in with above service account and leave message

請輸入暱稱 ( 最多顯示 6 個中文字元 )

請輸入標題 ( 最多顯示 9 個中文字元 )

請輸入內容 ( 最多 140 個中文字元 )

Please input verification code on left:

Cannot understand, change to another image

請輸入驗證碼