Edison.X. 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

以１０為底的對數，直接基於 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 報銷，一切從零開始。