[公告] 痞客邦新服務上線 每日星座運勢測算[公告] 痞客邦應用市集全新改版![公告] 痞客邦「應用市集」新 App 上架-iFontCloud Professional[公告] 痞客邦後台發表文章提供插入多張圖片新功能[公告]痞客邦新服務上線 部落客商店聚集就在《痞市集》

二年前瘋狂研究 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 留言(2) 引用(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

other options

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

請輸入驗證碼