前言

這是我曾幹過的傻事,若微積分、複變函式、線性代數不夠強的話,我希望別有人再重蹈我的覆轍除非本身工作是開發極為低階的數學函式,或本身之編譯環境裡,並沒提供低階的數學函式,本篇可做為一點參考。

本文將提到一些簡單的數學方式去應用低階的數學函式,這可以當作程式上之練習沒錯,我認為它也是個不錯的練習,但若對程式語言興趣缺缺,或覺得不實用,可跳過本文不予以理會。

非 CS 領域、對程式設計有一定熱衷程度時,通常會想要寫一些 ANSI C 裡面已經有的函式庫,比較簡單的可能像是 string.h (有些其實也不簡單了 ) ,甚至有機會碰到數值分析時,走火入魔時又發現,這些書都是基於現有的 math.h 去呼叫一些簡單的數學函式,如 exp、pow、Round function、三角函式、反三角函式、雙曲函式、複數係計算 等,這些太基本的函式,數值分析上的書都拿出來直接調用較多,於是走火入魔時,便有種異想天開的想法:怎麼做得比內建的 math.h 還快?

其實一些簡單的做法都是在考驗數學學得好不好,甚至有些數學題目是學校作業裡面就會提到的,或是有些在一些書上裡面也會提到的,但強調的是,他們的作法不是錯,而是效率可能不好,也可能有些全面性的問題沒考量到,以下針對一些常見的 math.h 函式做簡略說明。

 


 

Round function

 

fbas / abs 

fabs 之作法較沒問題,在 C 語言甚至覺得用 function-like macro 會較為方便,原因在於如此便不用區分為 abs 或 fabs 。

 

#define FABS(x)  ( (x)>0 ? (x) : -(x))

int abs(int x)
{
    return x>0   ? x : -x;
}

double fabs1(double x)
{
    return x>0.0 ? x : -x;
}

 

ceil / floor

至於 ceil 與 floor,一般人顯少討論到負數情況,連負數都考慮進去時大致是這樣

 

#define CEIL2(x) ((double)( (x)==(int)(x) || \
                  (x)<0.0 ? (int)(x) : (int)(x+1.0)))

double ceil2(double x)
{
    if( x==(int)x || (x<0.0) ) return (int)x;
    else return (int)(x+1);
}

 

#define FLOOR2(x) ((double)( (x)==(int)(x) || \
                   (x)>0.0 ? (int)(x) : (int)(x-1.0)))

double floor2(double x)
{
    if( x==(int)x || (x>0.0) ) return (int)x;
    else return (int)(x-1);
}

 

fmod

fmod 必須說明一下。一般而言,A / B,若 A, B 同號,沒這個問題;但若 A, B 為異號時,結果是 depends on machine,有些餘數會正,有些餘數會負,相對的,商之結果也會差 1。比較直接的做法,是將 A / B 取整數之商 C,再進行 return A - B*C 即可。但在不少情況,需要的是將商直接取 > 0,故建議分兩個函式寫會較佳。

 

double fmod1(double numerator, double denominator)
{
    return numerator - (int)(numerator / denominator) * denominator;
}

double fmod2(double numerator, double denominator)
{
    int quotient = (int)(numerator / denominator);
    numerator -= quotient * denominator;
    return numerator+ ( numerator<0.0 ? denominator : 0.0 );    
}

 

 

 


 

Exponential & Logarithmic Function

frexp / ldexp

有兩支函式較不建議寫,frexp / ldexp。frexp 是取得該浮點數之 Matissa 與 Exp 部份 (二進位表示),而 ldexp 則反過來,是給該浮點數之 Matissa 與 Exp,傳回表示之浮點數為何。不建議寫的原因在於,若這兩支寫下去,是否代表整台machine / compiler 之浮點數表示法都要自己搞、自己規劃?雖大多 comiler 都從 IEEE754,但一些為了表達更小、更大時,可能會有所不同。當然若是要了解 IEEE754 當練習,不反對,此處便不予以實做。

接下來所有的東西,全都由 Taylor series 發展出來。

exp

在寫 exp function 時,用 Taylor series 去展開可得如此

exponention.png  

這題目非常熟悉沒錯,但一般在實做上,不是決定要執行幾個回圈,而是先決定可接受最小誤差 eps 為多少,當上一迭代之解與下一迭代之解差額 <=eps 時,即可跳出回圈。我其實沒很喜歡用這種表達方式表示,exp, sin, cos 等用 taylor series 表達出來的,基本上都可以畫簡過,以 exp 為例,公式可化成下面這樣

exp(x) = 1 + x(1+ x/2! (1+ x/3!(... (x/n!)...)))

double expn2(double x, int n)
{
    double ret;    
    ret=1.0+x/n, --n;
    while(n) ret = 1.0+ret*x/n, --n;
    return ret;
}

 

這樣有幾個優點,運算速度略為快點 (是有一說,debug mode 會較慢,此處不予追究,有興趣請自行驗證),另精度「可能」比從前面算還精些。但這種寫法卻不非常實用,因無法估算出,exp(x) 必須在幾個迭代後,誤差才會小於 eps,若有估算大略迭代次數之方法,也請不吝提供。

該注意到的是,exp 將會是整個 math.h 效能好不好的瓶頸之一,用 Taylor series 方式展開,要達到 DBL_EPSILON 之精度,上述的寫法所費時間是 VC (別說 VC 了,gcc 可能都打不過了) 之 100 倍,即使將精度降低,也打不過 compiler 裡面寫好的 math.h。

modf

modf 是將一個浮點數,拆成整數及小數部份,若該浮點數為負數時,也是 depends on machine ,故也建議拆兩個出來寫,一個保證確定小數部份必大於零為佳。

 

double modf1(double x, double *intpart)
{
    *intpart = (int)x;
    return x-*intpart;
}

double modf2(double x, double *intpart)
{
    *intpart = (int)x;
    if(x-*intpart < 0) *intpart-=1.0;
    return x-*intpart;
}

 

log & log10

math.h 裡之 log 是以 e 為底之自然對數,即 log(x) 意指為 ln(x)。這部份也是關係到整個 math.h 效能好壞關鍵之一。目前我所知方法的確很糟。

nature_log0.png  

看到這公式應該傻眼了吧?積分?我能想到的,只有用 梯型積分法 、Simpson積分法、Romberg 積分法 去求 nature log,接下來什麼都不用玩了,最重要的兩個 math.h function - exp、log 都搞得亂七八糟的,這個 math.h 函式庫被我徹底搞爛了。

假設已做好了 log 之後,log10 也不是問題了,用到的是換底公式。

VCMPI000.png  

還想自己做 log2 也沒問題。


Power Function

pow

部份會發現,調用 pow 速度似乎會較慢,那是因為一般的練習題都是假設指數 (exp) 部份是正整數的情況,的確若指數部份為整數,自己寫會較佳。

 

double pow1(double base, int exp)
{
    double ret=1.0;
    if(exp>0){
        while(exp--) ret*=base;
    }    
    return ret;
}

 

但這也不是個好方法,有個 fastpower 的方式

 

double pow3(double base, int exp)
{
    double rst=1.0;
    while(exp>0){
        if(exp&1) rst*=base;
        base*=base,exp>>=1;
    }
    return rst;
}

 

上面這些都只考慮到,exp 大於零之情況,實際上 exp 仍可能小於 0 ,可再加上一點技巧做改善。

 

double pow3(double base, int exp)
{
    double rst=1.0;
    int abs_exp = (exp>0) ? (exp) : (-exp);
    while(abs_exp>0){
        if(abs_exp&1) rst*=base;
        base*=base,abs_exp>>=1;
    }
    return exp>0 ? rst : 1.0/rst; 
}

 

即使如此,仍沒討論到指數部份是浮點數的情況,how to ? 這裡用到了對數裡的換底公式。

VCMPI001.png  

化簡到這步之後,a 已知,ln(b) 上面用積分法做,exp (a* ln(b) ) 用 taylor 展開,速度慢到一個不行。不過我很好奇,為什麼學校老師都沒人出過這題目?這一題等於是要做 exp(a)、要做 ln(n) 積分法,是個練習的好題目。 

sqrt

根號的問題在網路上其實非常多的討論,我想源起是來自於 InvSqrt 問題,但我說現在 InvSqrt 真的已比不上 compiler 優化之速度及精度,大概也沒人會相信,存在的意義,我認為大概只剩寫硬體,優化能力不佳之 compiler 。

一般而言,sqrt 有直式開根號、牛頓法、巴比倫法去做,直式開根號不很建議去做,因其中會關係到搜尋演算法之問題 (你要猜數字是多少),而牛頓法化簡後,事實上就是巴比倫法之公式,只是巴比倫法開根號這方法比牛頓法早出來而已。

 

// Babylonian method, newton method
double sqrt1(double x)
{
    double x0 = 0.5 * x;
    double x1 = 0.5 * (x0 + x/x0);
    if(x==0.0) return 0;

    while( (x0 > x1+EPS) || (x1 > x0 + EPS) )
        x0=x1, x1 = 0.5 * (x0 + x/x0);
    return x1;
}

 


 

Other Functions

這篇文章,最重要的部份僅於此。看完一定會問:那三角函式、反三角函式、雙曲函式怎辦?說「僅於此」,是因為我很不想把一堆公式 po 上來,這樣倒不如回去再念念微積分更加實際,不是嗎?目前三角函式、反三角函式、雙曲函式之求法,「公開的」大多還是以 Taylor Series 為主 ,這部份到最後只會貼一堆 Taylor Series 之展開式,這部份到 wiki 上可以找到比我還詳細的資料,由我說實在沒什麼意義。

然而的確是有些「小技巧」可再磨過,如三角函式之奇偶性質、和差化積公式等。因一般而言,Taylor Series 在算三角函數時,必須額外注意「值域」之問題,這部份每個人習慣不同,有些人習慣調為 0~pi/4 處理,有些人習慣調為 0~pi/2 處理,也有人直接調成 0~2pi 處理,這裡不再贅述這些問題。 


How to Speed up?

回到一開始那年幼無知的夢想,怎麼做得比 math.h 裡面的東西還快?

一般在網路上,能夠 google 到的,大多都是用精確度換取時間,意思是他們程式碼本身之精度並不高!好笑的是,即使這些方法犧牲了很多精度,但仍打不過一些 math library ,不論是時間或精度都輸。有興趣的話可參考 dev-master 裡面的 一份討論,這份討論很精彩,值得細看。

關鍵在於,網路上大多之方法,依然跑不掉三個部份在組合 :Taylor Series 展開式、牛頓法少次迭代 (精度就變差)、以 bitwise hacker 操作 IEEE754 ,在 bitwise hacker 操作 IEEE754 方法中,常會出現一些讓人難以理解的 magic number ;有些是原作自己推出來,有些則是用電腦硬跑下去,期得到較小之誤差。

對於三角函式實做有興趣的話,可參考  這份網頁  ,裡面已把一堆 math.h 用到的函式包好,但也再提醒,會有誤差。若只是想找一般數學函式加速的話,可參考 OpenGL-Mathematics 這份網頁。目前沒實際用過,但據說效果還算不錯。

說了半天還是沒說到重點,怎麼做得比 compiler 內建的 math library 還快?

精度比別人差,速度能達到「相仿」就已非常勉強了,更不用到精度一樣,速度要比別人快。一些寫 compiler math library 的人,掌握了 compiler 內部函式之秘密,這點可能也不是一般人會知道的。何況,我想他們可能也不是用 Taylor Series 方式去求後面這些三角函式、反三角函式等,也可能是用 Fourier Transform 去求。

要比他們快,數學不好、沒念過一堆 paper,我想是不可能達到的,一個精度有限的 InvSqrt 就可以發一篇論文,想想 math.h 裡面有幾個函式,可以發幾篇論文?一篇論文最快算一年,那又要寫幾年?還是先去挖別人寫好的論文出來念念。但,念了可能還是寫不過 math library,因那些人很可能早念過這些論文。

真的要寫出無敵的 math.h 嗎?

這裡很不負責的結論是,別傻了。三角函式、反三角函式、雙曲函式、exp、log,能寫出一個與 math.h 相仿,我認為便可準備發一篇論文了。

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