前言

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

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

非 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 發表在 痞客邦 PIXNET 留言(10) 人氣()


留言列表 (10)

發表留言
  • Chen-Pang He
  • exp: 函數直接用泰勒級數跑比較慢,用 Pade approximant 來逼比較快。

    log: 其實是算 log(x + 1)。的確是用積的,不過要偷用 Legendre polynomials.

    矩陣函數當然沒有內建的,所以要自己寫。裡面的 maxnorm 相當於實數運算的最大絕對值。f(0)是精確值,所以絕對值越大越不準。絕對值太大就要把實數除以 2 (for exp) 或開根號 (for log)。

    裡面還有小弟最近才寫出來的 pow,請笑納!:)

    Eigen 矩陣函數的 C++ 源碼: https://bitbucket.org/eigen/eigen/src/97e16488639a/unsupported/Eigen/src/MatrixFunctions
    Eigen 官網: http://eigen.tuxfamily.org/
    GMP 官網: http://gmplib.org/
  • 謝謝您細心閱讀我的文章與回覆,進而請教幾個問題。

    (1) exp : 其實只是先再補充。 exp 上面我寫的是最初稿的,問題簡單時,上在求 exp(x) 可將 x 分解成整數/小數部份,整數部份可直接用 O(logn) 之 fast power,小數部份 exp(x) 才以 taylor series 在 |x| < 1 之情況下收斂性我實作已算 "不錯"。而您所提到 "Pade approximant" 我卻苦無合適資料下手,不知是否有何建議?

    (2) log(x) : 這個最後我不是用積分,其實會先化成 x = a * 2^b,其中 0.5 <= a < 1,再套入 ln(z) = (z-1) - (z-1)^2/2 + ....,想請問 Legendre polynomial 能否提點一二?

    (3) 您所提到的 Eigen 所附的三個 website,其中有兩份還蠻有名的,這部份請教是出自您手嗎??

    (4) 小弟文中所提的 pow 函式,指的是 double pow(double, double),也就是 2.34^4.83 這種數值,細閱您的回文,指的應是 matrix 裡之 power operation,在此先澄清誤解。

    最後謝謝您不吝提出之意見,感謝 :)

    edisonx 於 2012/09/02 01:08 回覆

  • Chen-Pang He
  • (1)
    Taylor series 是用多項式來逼近無理函數。Pade approximant 是用有理函數(兩多項式的商)來逼近無理函數。如果想要算 [m/n] 次也就是分子 m 次多項式,分母 n 次多項式的係數,要先算好 [m+n] 次 Taylor series 再用 Pade 去逼 Taylor。係數都是有理數,我們可以先算好精確的係數再 code 進去。

    Pade approximant 通常會比 Taylor series 快。2n 次的 Taylor series 其實就是 [2n/1] 次的 Pade approximant,不過 [n/n] 通常比 [2n/1] 準很多。
    (準就是快,因為我們看到誤差小於 epsilon 就/才停)

    (2)
    嗯嗯,用多項式逼近 log 有更快的算法噢!XD
    ln(1+x) = x - x^2/2 + x^3/3 - x^4/4 + ...
    ln(1-x) = -x - x^2/2 - x^3/3 - x^4/4 - ...
    ln((1-x)/(1+x)) = -2 (x + x^3/3 + x^5/5 + ...)
    y := (1-x)/(1+x) 則 x = (1-y)/(1+y) 即它是自己的反函數 :)

    不過,這只是多項式。化成有理函數會逼近得更快。
    log(1+x) = c[1]*x/(1 + c[2]*x/(1+ c[3]*x/(1 + ... + c[2m-2]*x/(1 + c[2m-1]*x/(1 + c[2m]*x)))))

    c[1] = 1
    c[2j] = j / (2 (2j - 1))
    c[2j+1] = j / (2 (2j+1))

    Gauss-Legendre quadrature 的根本(邏輯基礎)跟上面的 Pade approximant 一樣,不過實作差很多。Gauss-Legendre 是梯形積分,不過每個梯形的高不一樣長就是了。

    [維基百科條目,上英下中]
    http://en.wikipedia.org/wiki/Gaussian_quadrature
    http://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%B1%82%E7%A7%AF

    理論上它們收斂的速度一樣,不過對矩陣而言是 Gauss-Legendre quadrature 比較準啦。對實數的話,您可以用 wxMaxima 之類的東西跑跑看。Sorry, 小弟最近剛開學比較沒空。(謎:醫二不是很閒嗎?)

    (3)
    我們目前主要參考 Nicolas J. Higham 的論文。指/對數是 Jitse Niesen 寫的,我給它加上對 long double 的支援。

    矩陣.pow(實數) 是我最近寫的,目前有一點效能不太好的地方就是:如果要算同一個矩陣的次方很多遍,像是要算 A.pow(0.7),也要 A.pow(1.1),會重做一遍 Schur 分解 (A = U T U^H, U 是么正(unitary)矩陣,T 是上三角矩陣),因為三角矩陣才有快速的 T.pow(小數) 可以用。我現在正在把它改得可以儲存分解結果,讓使用者可以直接把它 call 出來。值 25 次矩陣乘法的 Schur 分解就不要再勞駕它了!XD

    (4)
    沒有誤解啊,實/複數是 1x1 矩陣,當然吃太飽也可以把複數看成 2x2 實矩陣。:)
  • Chen-Pang He
  • 啊,對吼,還有實數指數(實數指的是指數噢!底數隨便,連矩陣都可以 XD)。複數指數就用 log then exp 吧,隨它去。

    pow(1-x, p) = 1 + c[1]*x/(1 + c[2]*x/(1+ c[3]*x/(1 + ... + c[2m-2]*x/(1 + c[2m-1]*x/(1 + c[2m]*x)))))

    c[1] = 1
    c[2j] = (-j+p) / (2 (2j - 1))
    c[2j+1] = (-j-p) / (2 (2j+1))

    言歸正傳,我還是覺得我們拼不過 math.h 啦!如果有 90% 的速度的話,或許我們就可以出論文了? 不過 pow 搞不好有機會可以贏 math.h 噢,您可以拿去測測看 XD,因為它搞不好是用 log then exp。我是沒多少信心啦,所以我們目前還是調用內建函數來跑矩陣的主對角線。

    GCC 在 Windows 下要有 vectorization 或 multi-threading 才有可能比 MSVC 快(吧?)。Eigen 在有 OpenMP 的情況下會用它來做 multi-threading, MinGW 有內建,不曉得 MSVC 裡面有沒有。小弟是窮光蛋/自由軟體擁護者,非必要不買 proprietary software;現在除了交作業有時候要用 Office 之外,幾乎都在 Fedora 下過活。XD
  • 感謝你的留言與細心指教,真的蠻多東西沒碰過,(可能我數值分析的書講太淺?) 謝謝不吝分析您的知識與看法。

    最後小提一下下,math.h 蠻多函式在套入公式前,會直接取 IEEE754 欄位(也就是 frexp) 進行切割後再做運算,加速收斂速度。

    edisonx 於 2012/09/14 20:05 回覆

  • Chen-Pang He
  • 嘿呀,像是 log 的話,機器可以先「背好」log(2) 的值,我們可以先算好再給他 literal 進去。其實算 log 的時候,我們是算 log(1+x)。不過 1+x 的範圍最好不要取(很顯然的) 0.5~1.0 或 1.0~2.0,因為 log(1.1) 比 log(0.55) 好算,log(0.9) 也比 log(1.8) 好算。成本最低的區間 c~2c 會在這兩者之間。如果 x>0,log(1+x) 會比 log(1-x) 準,所以 c 應該會比 2/3 大。我猜應該不會離 sqrt(2)/2 = 0.7071...太遠。

    算浮點數跟算矩陣最大的不同是:矩陣乘法很貴,加法跟讀取根本不算什麼。所以成本分析很簡單。浮點數在做成本分析的時候,加、減、乘的成本其實差不多,讀取的成本也不能忽略。這部份您可能要另請熟悉硬體的高明了,小弟習慣躲在 C/C++ standard 很久了。XD

    不過算浮點數比矩陣輕鬆的部份在於他被存在 mantissa * 2^exponent,所以可能不用像矩陣一樣擔心收斂性能。

    小弟寫出了還算可以的 exp,運行時間大概是內建函數的兩倍多一點。請笑納!

    #include <math.h>
    #include <stdio.h>
    #include <stdlib.h>
    #include <time.h>

    double pade1(double x)
    {
    return (2+x) / (2-x);
    }

    double pade2(double x)
    {
    double u = 6*x;
    double v = x*x + 12;
    return (v+u) / (v-u);
    }

    double pade3(double x)
    {
    double x2 = x * x;
    double u = x * (x2 + 60);
    double v = 12*x2 + 120;
    return (v+u) / (v-u);
    }

    double pade4(double x)
    {
    double x2 = x * x;
    double u = x * (20*x2 + 840);
    double v = (x2 + 180)*x2 + 1680;
    return (v+u) / (v-u);
    }

    double pade5(double x)
    {
    double x2 = x * x;
    double u = x * ((x2 + 420)*x2 + 15120);
    double v = (30*x2 + 3360)*x2 + 30240;
    return (v+u) / (v-u);
    }

    double pade6(double x)
    {
    double x2 = x * x;
    double u = x * ((42*x2 + 10080)*x2 + 332640);
    double v = ((x2 + 840)*x2 + 75600)*x2 + 665280;
    return (v+u) / (v-u);
    }

    double compute(double x)
    {
    double norm = abs(x);

    if (norm < 1.100349098426267e-5)
    return pade1(x);
    else if (norm < 2.401864567319500e-3)
    return pade2(x);
    else if (norm < 1.495585217958292e-2)
    return pade3(x);
    else if (norm < 1.122050472231552e-1)
    return pade4(x);
    else if (norm < 2.539398330063230e-1)
    return pade5(x);
    else /* (norm < 5.676056633477773e-1) */
    return pade6(x);
    }

    double expi(int i)
    {
    double res, tmp;

    if (i>=0) tmp = 2.718281828459045;
    else {i=-i; tmp = .3678794411714423;}

    res = (i&1) ? tmp : 1;

    while (i>>=1) {
    tmp *= tmp;
    if (i&1) res *= tmp;
    }
    return res;
    }

    double expd(double x)
    {
    int intexp = floor(x + 0.5);
    x -= intexp;
    return compute(x) * expi(intexp);
    }
  • Chen-Pang He
  • 啊啊,不小心 include 太多東西了
    #include <math.h> /* 有必要,因為要用 floor 跟 abs */

    剩下的是我自己測試時 include 進去的。=.=
    縮排全部不見了,可讀性瞬間降低。(汗)
  • Chen-Pang He
  • 不知為何 pade1 ~ pade5 的誤差都頗大,不如理論預期。所以建議把它們刪掉,直接把 pade6 的內容搬進 compute 就好。這樣大概慢個 5% 但是準很多。
  • 這份 code 確實寫得漂亮,Pade approximant 我可能需要找時間補一下 XD。 Pade 竟只要六次就保證收斂,這點的確讓我有點吃驚,最後感謝您的分享 :)

    edisonx 於 2012/09/16 16:34 回覆

  • Chen-Pang He
  • 唉呀,我剛 benchmark 了一下,發現內建的 exp 無論參數多大,算的速度都一樣。所以我懷疑他是用其他演算法。然後我猛然想到了:

    算矩陣因為精確度的問題,是直接取 exp 的 Pade。算實數好像用 coth 比較快,沒記錯的話是快一倍。(驚)

    我等一下去跑一下 coth 的係數跟容忍值,回頭見!:)
  • Chen-Pang He
  • 太有趣了,找到癥結了:要先估算結果是 2 的幾次方。然後因為自變數縮到正負 log(2) / 2 之間,所以 Pade 5 次就夠了。我這邊用 gcc 測,讓它跑 2^24 次:

    無: expd 1320ms exp 750ms
    -O: expd 800ms exp 650ms
    -O2: expd 750ms exp 190ms
    -O3: expd 740ms exp 180ms
    -Os: expd 720ms exp 650ms

    用 define 是為了等一下來玩 log 或 pow :)

    #include <math.h>

    #define C1 0.693359375 /* C1 + C2 should represent log 2 to */
    #define C2 -2.12194440054690583e-4 /* more than machine precision */
    #define LOGB2E 1.442695040888963

    double expd(double x)
    {
    const int b[] = { 30240, 15120, 3360, 420, 30, 1 };
    double u, v, z = x * LOGB2E;
    int n = (z < 0) ? z - 0.5 : z + 0.5;

    x = x - n * C1 - n * C2;
    z = x * x;

    u = ((b[5]*z + b[3])*z + b[1]) * x;
    v = (b[4]*z + b[2])*z + b[0];

    return ldexp((v + u) / (v - u), n);
    }
  • 這段程式碼真是讓人吃驚...

    edisonx 於 2012/09/16 23:18 回覆

  • Chen-Pang He
  • 把 exp 和 log 寫在一起嘍,因為都會用到 C1 跟 C2

    #include <math.h>

    #define C1 0.693359375 /* C1 + C2 should represent log 2 to */
    #define C2 -2.12194440054690583e-4 /* more than machine precision */

    #ifndef M_LOG2E
    #define M_LOG2E 1.442695040888963
    #endif

    #ifndef M_SQRT1_2
    #define M_SQRT1_2 0.7071067811865476
    #endif

    double expd(double x)
    {
    const int b[] = { 30240, 15120, 3360, 420, 30, 1 };
    double u, v;
    int n;

    v = x * M_LOG2E;
    n = (v < 0) ? v - 0.5 : v + 0.5;

    x = x - n * C1 - n * C2;
    v = x * x;

    u = ((b[5]*v + b[3])*v + b[1]) * x;
    v = (b[4]*v + b[2])*v + b[0];

    return ldexp((v + u) / (v - u), n);
    }

    double logd(double x)
    {
    const double nodes[] = { 0.0254460438286207, 0.1292344072003027,
    0.2970774243113014, 0.5000000000000000, 0.7029225756886985,
    0.8707655927996972, 0.9745539561713792 };
    const double weights[] = { 0.0647424830844348, 0.1398526957446383,
    0.1909150252525594, 0.2089795918367346, 0.1909150252525594,
    0.1398526957446383, 0.0647424830844348 };
    double r, t;
    int i, n;

    t = frexp(x, &n);
    r = 0;

    if (t > M_SQRT1_2) {
    x = t - 1;
    }
    else {
    --n;
    x = t * 2 - 1;
    }

    for (i = 0; i < 7; ++i) {
    r += weights[i] * x / (1 + nodes[i] * x);
    }
    return n * C2 + r + n * C1;
    }
  • 請問 nodes[] , weights[] 這些常數和公式是 k 過 Pade 就該懂的嗎?是的話我想我該再 k 一次 Orz。

    edisonx 於 2012/10/01 02:23 回覆

  • Chen-Pang He
  • 噢噢,nodes[] 跟 weights[] 是拿來做積分的。
    http://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%B1%82%E7%A7%AF

    程式碼裡面的值有位移過:nodes 存的是 (x+1)/2,weights 存的是 w/2。這樣會比原公式準一點點。

    這兩個的演算法都直接拿 Eigen 算(三角)矩陣的來用。:)

您尚未登入,將以訪客身份留言。亦可以上方服務帳號登入留言

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

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

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

請輸入左方認證碼:

看不懂,換張圖

請輸入驗證碼