二年前瘋狂研究 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 報銷,一切從零開始。

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