天天看點

第七章 數學

說明:本文翻譯自《TangoRefMan_Sep_1_2008》

由于本人是程式設計初學者,對很多程式設計概念不是非常熟悉,程式設計經驗不多,再加上英語水準不高,翻譯純屬一個D語言愛好者實驗之作,很多錯誤在所難免,還請讀者見諒。另外,如果你發現本文有不當和錯誤之處,還請多提寶貴意見。  

第七章 數學(Doing the Math)

介紹(Introduction)

計算機最初的構想是作為進行數學運算的裝置。早期的電腦花費大量時間用于求解方程。盡管現在工程和科學計算隻是計算世界的一小部分,然而計算技術的發展留下來一個寶貴的财富:那就是現在幾乎所有的電腦硬體具有超高的計算準确性和非常快的計算速度。遺憾的是,大多數程式設計語言很難讓程式員充分利用硬體的這些特性。一個顯著的問題是缺乏文檔;甚至對許多數學程式員來說,浮點運算仍然是神秘的事物。

在大多數的程式設計語言中,對數學地支援僅限于單個檔案。例如,在C和c++,所有的數學函數都包含在“math.h”中。考慮數學領域的範圍,以及數學與計算科學的曆史聯系,這個有限的存在是幾乎是一種侮辱。這種安排也意味著深奧的功能與簡單的和廣泛應用的功能混雜起來,使得該領域的新人更加感到撲朔迷離。在tango數學庫中,基本的“高等學校”數學函數從低級數值分析函數和更少熟悉的特殊函數中被分離出來。

IEEE浮點運算(IEEE floating-point arithmetic)

tango的數學庫更多地涉及浮點運算。要了解它,必須對IEEE浮點算法有一個基本的了解。從根本上說,它是一個無窮多的實數在一部分位元組上的映射。作為一個32位浮點數,僅有40億個不同的數字可以被正确表示。正是在這樣一個讓人沮喪的狹小的表示空間中,IEEE浮點算法做了一個不同尋常的維護幻想的好工作,使數學實數可以正常工作,不過重要的是要明白什麼時候打破了幻想。

第一個缺點來自這些可表示數字的分布。IEEE數軸與數學上的數軸迥然不同。

+  +-----------+------------+.. +.. +----------+----------+ +     #

-inf -real.max  -1    -real.min   0   real.min   1   real.max inf  NAN

注意到一半數量IEEE數軸在-1和+ 1之間。當使用80位實數,有超過0.5到1之間一千倍的從0到0.5間的數。(這不是一個純粹的對數關系,那就是:在1.5和2之間存在與1和1.5之間同樣多的可正确表示的數)。這對于精度具有重要意義:有效的精度令人難以相信地接近零。利用這一點在-1和+1間的數看作分開的會讓一些例子被覆寫。

同時也要注意這些特殊的數:+ / -無窮大,所謂的“denormals”或"subnormals"在+/-real.min 和 0之間,以減少的精度來表示,最後神秘而令人難以置信的有用的“NaN”(“Not-a-Number”)将在下面讨論。

同時也要注意這些特殊的數:+ / -無窮大,所謂的“denormals”或"subnormals"在+/-real.min 和 0之間,以減少的精度來表示,最後神秘而令人難以置信的有用的“NaN”(“Not-a-Number”)将在下面讨論。

 

 第二個重要概念是需要注意取整誤差(round-off error)。由于大多數數學上的數字不能完美表示,隻好用最接近的可描述數字來代替。是以,幾乎所有浮點運算會遭受一些取整(舍入)誤差。不幸的是,取整誤差會在漫長的計算中積累。在Tango數學庫執行浮點計算的一些函數采用了不同尋常的方式,以最小化取整誤差。基于同樣的原因,幾乎所有函數均隻有real的精度。

  

有必要知道哪些運算可以避免取整誤差。任何可以存儲在int型中的數可以在double中精确表示,任何可以存儲在long型中的數可以在real中精确表示。任何與int進行運算時精确的算法采用浮點時也是精确的(如,27.0/3.0 = = 9,精确)。與任意的2的幂進行乘、除也是精确的。

了解這些對于進行相等測試時具有重要意義:不要使用==進行精确的相等測試。該測試通常隻能得到近似相等,它帶有一些舍入誤差的修正值。Tango庫中包括一個函數feqrel()可用于該目的的,将在下面讨論。

錯誤處理

當一個提供的參數超出函數可表示的領域時,函數會怎樣做呢?例如一個最簡單的情況sqrt(-1)。有如下一些可能性,各有優缺點:

使用一個輸入的契約。建立調試版本時一個斷言被産生,進而幫助人們快速發現程式錯誤。不幸的是,這種行為在建立發行版本時是未定義的,除非另一個解決方案被很好的使用。另一個問題是,有效的定義域會變得非常複雜,它會不切實際的期望使用者代碼避免所有的違反契約情況。

抛出一個異常。這會導緻代碼膨脹,并強迫使用者把代碼弄亂放入try..catch語句塊中。

傳回一個NaN。一個Not-a-Number意味着它不可能提供一種明智的答案;例如,當要求CPU計算0.0/0.0,它會創造一個NaN。NaNs被硬體非常好地支援,并且調試時非常有用,不過會被大多數程式設計語言忽略。

在Tango.math中,所有的函數在表示無效輸入時均傳回NaN。此外,契約僅用于在明确知道有一個程式的邏輯錯誤時使用(例如,對于統計函數“statistical functions”來說,所有可能性必須在0 . .1之間)。不過它會盡可能的做得比這更好一些。IEEE NaN可以包含一個有效負載——一個整數存儲在NaN内,該整數可以用來解釋NaN來自哪裡和它為什麼會産生。D也許是第一個使用這個性能的流行程式設計語言,是以必須制定一個ABI。如果可能的話,在可能的地方,tango庫NaN包含一個有效負載解釋它在哪裡被産生。這個有效載荷被I/O函數顯示,也可以被庫函數取回。注意,在許多處理器上包含NaN的計算通常都比較慢,是以更可取的是要限制他們使用到異常情況下。

math.Math

在math.Math中,最常見和通用的處理浮點數的數學函數都可以找到。幾乎所有的函數都會被學過高等數學的人所熟悉。這些函數可分為三類:乘方、三角、舍入。

(常量)Constants

與同等的數學常量有一些微妙的差異。注意,PI不是數學數字pi;它隻是在可表示的範圍中最接近PI的那個數。是以,例如,sin(pi)= 0,但罪(PI)= 2e-20。這個取整誤差的原因不在正弦函數sine上,而在于PI的表示精度。(實際上,情況更複雜的:在x86平台,sin(x)在abs(x) > pi/2的無論何時都會傳回不正确的值)。

乘方、根和絕對值(Powers, roots, and absolute value)

abs: 可以是integer, real, complex,和 imaginary 參數

sqrt: real 或 complex,。 real版本決不會傳回一個 complex 結果。

cbrt :立方根(cube root)。 僅有real 型。

exp, log, expm1, log1p.

exp,log是指數和自然對數函數.對于大的正數和負數exp()是精确的,不過當任意參數非常接近0時當然會傳回1.0。像先前讨論的一樣,IEEE數軸表示接近0的數要好于表示接近1的數。expm1(x)等價于exp(x)- 1.0,但它儲存精度幾乎降為零。同樣的,log1p(x)是一個相當于log(x)+ 1當x的值接近零時的值。在适當的時候,用這兩個函數代替exp()和log()是一種提高精度的簡單的方法。當然,x是一個任意的大負數x時,expm1(x)恒為- 1。

pow -- 乘方, 高速計算整數的乘方

hypot -- 計算直角三角形斜邊的長度,應用畢達哥拉斯規則(Pythagoras' rule),但用了一個避免溢出的方法。

poly -- 計算矩陣多項式(polynomial )。

conj -- 共轭複數的計算 。

max, min, maxNum, minNum

最大值和最小值是出乎意料的複雜。如果一個參數是一個NaN将會發生什麼情況呢?根據754R标準,如果一個參數是NaN而另一個是正常數,該數應當被傳回(對于檢測一個函數的範圍來說,該行為是有用的)。是以,這些函數被命名為minNum 和maxNum,以提醒不論何時應盡可能地傳回一個數。遺憾的是,它的速度要比樸素的執行要慢一點。max and min 是 maximum和 minimum 函數工作在任意類型 時的版本;若任何參數是一個NaN,則傳回結果是未定義的。

三角函數

對于實數(real)、複數( complex)和虛數(imaginary )參數,全部用弧度(而不是度)。

cos, sin, tan

acos, asin, atan

cosh, sinh, tanh

acosh, asinh, atanh

即使對于那些不熟悉的複數的人,也值得知道下面的聯立的(同時的)快速計算正弦和餘弦的方法。

creal z = expi(real x);

real c = z.re;    // = cos(x)

real s = z.im;    // = sin(x)

在大多數處理器(包括x86),這是一種固有的操作,速度是分别計算正弦和餘弦的兩倍。

注意就像 exp(x), cos(x) ,當x很小時通常傳回1。cos(x) - 1形式的表達式将會用 cos(x) - 1 = -2 * sin(x/2)2 形的結果重寫。

sinPi, cosPi, atanPi

要精确地利用度代替弧度進行三角函數運算,可以使用sinPi(x / 180)。自從360可以被浮點數正确地表示,但pi不能,這些函數仍然對大的轉數保持準确。(例如,sinPi(1e30)正确地傳回0,但sin(1e30 * PI)的傳回值是無意義的。

函數sec = 1/cos, cosec = 1/sin、cot = 1/tan,它們曾經在三角學中是重要的,但在對數被介紹之後就變得不那麼重要了(http://mathworld.wolfram.com/Secant.html)。他們并不包括在tango.math中,因為它們用手工計算是微不足道的,還因為sec()這個名字和單詞seconds(秒)非常相似,會造成不必要的混亂。

atan2(y, x)

這是用來轉換直角坐标(x,y)成極坐标(angle, radius(角,半徑)),這對于複數運算是基本操作。(它也是适合用在零的符号特别重要的場合上(本句不夠準确))。

舍入到整數(Rounding to integers)

将一個實數轉換到最接近的整數是非常普遍的需要。有一個小小的困難:當存在一個難以取舍的平局情況怎麼辦呢?3.5應當舍入到3還是4 ?為了相容C語言,cast(int)通常向上舍入。不幸的是,在許多處理器上這個操作出乎意料的慢,因為它需要CPU的的舍入的方式被改變,随後又重新恢複。(它是如此緩慢以緻于Pentium 4專門包括一個特殊的優化器以緩解速度的懲罰,SSE4指令集進行進一步的優化)。然而在大多數應用程式中,采用何種舍入方式不是很重要的。函數rndint和rndlong是cast(int)和 cast(long)的高速替代者。它們的決定權(tie-breaking behaviour)取決于舍入模式,稍後将進行讨論。

還提供了一些另外有用的從實數到整數轉換的模式

floor

ceil

round

trunc

值得注意的是這些函數傳回實數(real),而不是整數( integers)。要儲存它們為 integral類型,仍然還要使用rndint或rndlong。

math.IEEE

math.IEEE是低級函數的家,它與計算機硬體的密切關系要超過與數學的關系。要讓這些函數變得有用,需要對IEEE算法有一定的了解,不過傳統上它們通常被給予非常晦澀的名字,在數學庫之外很少被使用。人們希望他們會更廣泛用于D,因為他們既能大大提高計算精度和健壯性,他還可讓一些有趣的技巧成為可能。

檢查近似的相等(Checking for approximate equality)

操作符==在大多數相等測試中通常隻進行精确地相等比較,因為哪怕是一個微小的取整誤差就會導緻測試失敗。相反的,一些足夠接近(close enough)形式的測試必須被使用。不幸的是,很容易就能讓寫出的這個測試不正确。例如,很明顯abs(x - y)<=0.000001是不正确的,比如,它意味着1e-10近似于7.8e-300,但這并不意味着1.000000000000000000000001e60近似等于1e60。取而代之的是,我們需要一個在IEEE數軸附近的測試而不是“在數學實數數軸附近”的測試。

分析師通常所說的對待取整錯誤的态度—"units in the last place" (ulp)(把誤差留在最後的地方),但這僅僅是誤差很小的時候才有用。一個更直覺的措施來統計舍入錯誤的位數(或做另一個選擇,保留準确的位數)。雖然這聽起來有點複雜,它可以充分地利用二進制表示浮點數的優勢來有效地執行。在編譯時指定允許的最大舍入位數,可以使測試速度非常快——隻是稍稍慢于正常的==比較。為确定清爽的舍入誤差是多大,使用tango.math.IEEE中的準确函數feqrel() 。

feqrel ( D特有)

Returns the number of significand bits which are equal in both x and y. Formally, real.mant_dig - ceil(log2) of the number of times nextup() would have to be called to move from the smaller number to the larger one.

IEEE 異常狀态标志(The IEEE Exception Status Flags)

  

所有IEEE-compiliant處理器包含特殊狀态位,用來訓示程式可能想了解的當“異常(weird)”的事情發生時狀态。例如,ieeeFlags.divideByZero(IEEE标志.被零除)告訴我們一個數被零除時一個無窮大發生了。在計算的末尾僅檢查一次,就可以避免相比成千上萬次的失敗比較。

  

這是一個可以探測到的異常狀态的清單:

無效(invalid):設定為當任意NaN形成時。這個狀态可能發生于∞-∞,∞×0,∞×0,0÷0,∞÷∞,∞%∞,x%0(x是任意數),以上這些情況幾乎總是預示着一種程式設計錯誤。

被零除(divisionByZero): 設定為被零除正負無窮大發生時,這通常預示着一個程式設計錯誤,但不是經常發生; 某些類型的計算當被零除情況發生時傳回正确的結果。 (如, 1/(1+ 1/x) == 0當 x == 0時)。

溢出(overflow):設定為把兩個非常大的數字相加或相乘,它們的結果大于real.max時。

下溢(underflow):當把兩個非常微小的數相減或相除時,它們的結果接近零時就産生下溢。下溢幾乎從不産生問題。

不精确(inexact):幾乎所有浮點運算操作都設定這個标記!對精确性問題,已經進行了讨論。

IEEE取整模式(The IEEE Rounding Modes)

取整被函數setIeeeRounding和getIeeeRounding控制 。可以設定四種取整模式。預設模式舍入到最接近的數,在統計學上是最準确的,但缺少直覺。在平局的情況下,舍入到偶數。

取整模式             rndint(4.5)             rndint(5.5)               rndint(-4.5)   

舍入到最接近的數  4                              6                          -4                 平局情況下舍入到偶數

向下舍入               4                             5                          -5   

向上舍入                5                             6                         -4   

向0方向舍入           4                             5                         -4   

通常有兩種類型的原因要改變取整模式。一是簡單的檢查數值穩定性:當取整模式被改變時如果計算結果産生難以控制的不同結果,它(譯者按:它指改變取整模式)可以作為一個清除符号可以免受取整誤差帶來的麻煩。第二個原因是為了實作區間四則運算。但最明顯的用途是暫時控制rndint的意義,這樣它才能準确複制在一個内循環中的cast(int)的工作情況。

操作指數和有效數字的函數(Functions for manipulating exponent and Significant)

由于IEEE采用了基于2的指數,用2乘或除一個數是一種精确的操作(完全不會産生任何舍入錯誤),當然,這個結果也并不會産生上溢或下溢。利用這一事實,計算的精度有時可以得到提高,或者延展其使用範圍。函數ldexp()和frexp()在計算結束時,允許從一個數或它們的恢複中去除所有的因素。

frexp -- 擷取指數和有效位數。

ldexp -- N * 2y

ilogb, logb

scalbn -- x * 2n

nextUp, nextDown, nextafter

nextUp(x)和nextDown(x)傳回 IEEE數軸上相鄰的點。他們已經加入了754R IEEE标準草案。nextafter()的使用是沮喪的,許多IEEE标準委員會成員現在認為這是個錯誤。

ieeeMean ( D特有)

我們在前面看到從IEEE數軸上存在的0和1間的數是1和2間存在的數的幾百或幾千倍。試圖用二分法(binary chop)獲得對數性能的分而克服(divide-and-conquer)算法是有問題的。考慮查找f(x) = 0中的x的值的問題,使用平分法把區間一分為二。如衆所周知的0 < x < 2,一個樸素的二等分會把該區間分成[0,1]和[1,2]。糟糕的是這是錯誤的,因為在區間[0,1]上包含了原區間上超過99%的可表示數字! ieeeMean(x,y)給出了表示範圍的中點,在所有的情況下提供良好的性能。注意:這個函數是沒有被IEEE标準描述。

signbit

除了要差別正負零外,與x > 0相同。

copysign

可以提高精度和速度的輔助函數(Helper functions that can improve accuracy or speed)

fabs --内在的。被 math.Math.abs使用。

  fma --代替“融合乘加”("fused multiply add")。(x * y)+ z。有些cpu(但不是x87)可以在一條指令中執行,這增加了準确性,因為它僅包含一個舍入(取整)誤差。

  expi—sin + cos.。在x86機器上這是一條指令。被exp(ireal)使用。

處理特殊實數類型的函數(Functions for dealing with special types of reals)

isNormal

isSubnormal

isInfinity

isNaN

  這些函數存在的主要為了相容c。注意isNaN(x)等價于x!<>=0,除了NaNs信号産生時不抛出異常外,isInfinity(x)與abs(x)= = x.infinity一樣。

  isIdentical(D特有,判斷同一性)

  除了正負零被看作不同外,相同作為相等。如果他們有相同的有效載荷,NaNs被視為相同。

NaN有效載荷(NaN Payloads)

根據IEEE 754标準,有效載荷(payload)可以存儲在NaN的尾部。該有效載荷會包含 NaN怎樣和為什麼産生的資訊。從曆史上看,幾乎沒有任何程式設計語言使用過這個潛在的強大功能。在tango中,該有效載荷由一個正整數組成。有一個困難是每當轉換到一個更小類型的情況發生時,某些位會從有效載荷中丢失。IEEE标準沒有明确說明此時發生了什麼事。

NaN — 帶“有效載荷 ”創造一個NaN,該“有效載荷”是一個ulong型數字。

ulong getNaNPayload(real x)——傳回整數載荷。注:如果縮小轉換已經發生,高序位可能已經發生了變化。決不要存儲一個指針作為一個整數載荷在一個NaN中。垃圾收集器将無法找到它!

 math.BigInt

 

對于可靠的應用程式來說,内建的整數類型精度不足,如某些應用高精度的數學計算和密碼程式等。tango使用BigInt類型提供任意精度的算法,BigInt能夠執行數千位精度的整數計算。一般來說,BigInt可以作為int和long型的插入代替者。要做到這一點,BigInt使用數值語義理論為指導,在D1.0 中用寫時複制(Copy-on-Write)實作。不幸的是,這個結果會産生很多臨時堆對象(在D2.0中不會發生這個問題)。

BigInt旨在提供優良的性能,用小一些的數字來緩和數位(勝任幾千個十進制數位)。高度優化彙程式設計式是用于最廣泛應用的處理器。對極其巨大的資料,或較為不通用的處理器,可以考慮用GMP(http://www.gmplib.org/)。

math.Bracket(支架、括号)

幾個重要的數值算法都包含在Tango中。數學。Math.Bracket包含查找根和求單一參數實數函數的極值。他們被設計用于較好地解決高等非線性函數。

FindRoot

經典的一進制根查找問題就是:給定a,b和一個實數函數f(real),計算f(x)==0(a<=x<=b)中的x。如果您可以提供一個範圍[a..b],這樣f(a) 和 f(b)有相反的符号,則要求的根可以在有效的範圍之内找到。(“有效”是指以盡可能少的調用f())。

Math.Bracket.findRoot使用一個原本基于TOMS478的算法(Alefeld et. al.)。Math.Bracket.findRoot使用一個原本基于TOMS478的算法(Alefeld最初et.。),盡可能用三次插值方法,在必要的時候退化到較慢收斂方法。顯著的速度和穩健性增強了。TOMS478與幾乎所有其他公衆可用算法遭受極其差的糟糕性能(“行業标準”算法由于P.R.布倫特可能要求為80位實約 40000次調用f());該實作避免了這樣的問題,使它比平均速度快1.5倍,在比最糟的情況快200倍。典型的10 - 15次調用f()需求達到全部80位精度。

例如,計算x5 = 0.2中的x,已知x在0到5之間,使用:

real root = findRoot( (real x) { return pow(x, 5) - 0.2; }, 0.0, 5.0);

如果該算法的支架(bracket)沒有包括零會工作得更快。通過再次觀察IEEE數軸,你可以看到,如果你可以指定一個範圍1e-5 ..5.0,你将給該算法比起你給一個包括0的範圍少做很多的工作。

注意不是所有的查找根據的問題都能夠用這種方法解決;要提供 f(x)上相反符号的兩點可能會變得複雜(有時f(x) 僅接觸x軸)。

FindMinimum

 

在一個次元的最小值(Minimisation)普遍使用布倫特(Brent)的算法,該算法結合了抛物插值法和平分法。

  在一個次元的最小值(Minimisation)普遍使用布倫特(Brent)的算法,該算法結合了抛物插值法和平分法。如 J. Crenshaw (http://www.embedded.com/story/OEG20010628S0046)說的,布倫特的算法并不是最佳的,将來tango的發行版中可能會有一個結構來取代它。這個函數僅僅找出機器精度的一半的最小值,因為取整誤差會讓更高的精度毫無意義(雖然不是經常,例如abs(x)有一個非常精确的确定最小值)。

  

如果你需要找到一個函數f(x)的最大值,隻要找到-f(x)的最小值。

注意,提供兩個點a和b給minimum做為托架,該算法也需要一個第三點c在[a..b]之間并滿足f(c)<= f(a) 和 f(c)<=f(b)。依賴于你的函數,查找三個這樣的點會極其複雜。

math.Probability(可能性、機率)

統計分布函數被廣泛地應用于數理統計理論,經濟學,生物學,以及整個實驗科學與工程學。Tango提供高精度最重要的機率分布工具。我們之前已注意到,介于0到0.5之間存在着0.5到1之間的數千倍的可表示的數。當表示極端的可能性時這可能會是個問題,因為通常存在同樣多的幾乎必然發生的事和幾乎一定不發生的事。要保留精度,函數被提供使用補助機率q,定義為q = 1 - p;對于接近1的機率使用蒙特卡羅版本(Compl)的分布函數。(這完全類似exp() 和expm1()之間的關系)。

在Tango中包含下面的累積分布函數。在每種情況下,補充(complement,字尾Compl),相反的(inverse,逆,字尾Inv)和補充和逆(complement inverse,字尾ComplInv)被提供。

MATHEMATICAL  TANGO NAME

 normalDistribution

Student's t  studentsTDistribution

f  fDistribution

chi2 chiSqrDistribution

 poissonDistribution

 binomialDistribution

 negativeBinomialDistribution

 betaDistribution

 gammaDistribution

  逆函數是用來統計資料構造置信區間(confidence intervals,也可以是可靠區間)。

  

  幾乎所有的這些函數的實作都是通過伽瑪(gamma)和β( beta )特殊函數,在下面将進行讨論

  數學的特殊函數(Mathematical Special Functions)

  數學術語的特殊函數(Special Functions)涵蓋大範圍的超越一般常識的函數。這是罕見的同時用在一個單獨的應用程式中查找倍數特殊函數(find multiple special functions)。

math.ErrorFunction(錯誤函數)

錯誤函數和正态分布函數的從根本上來說是相同的函數,從比例因子中分離出來。

erf, erfc

錯誤函數和互補錯誤函數。它們出現在最常見正态分布的內建中,也被稱為“鐘形曲線”(bell curve)或“高斯分布”(Gaussian distribution,盡管高斯并沒有發明它!)。錯誤函數和互補錯誤函數。它們出現在最常見正态分布的內建中,也被稱為“鐘形曲線”(bell curve)或“高斯分布”(Gaussian distribution,盡管高斯并沒有發明它!)。erf()接近0時有相對高的精度,而erfc對大正數有相對高的精度。(這類似于前面提到的expm1() 和 exp()的關系)。利用下面的結果,它給任意x盡可能獲得相對高的精度:

erf(x) = 1 - erfc(x)

erf(x) = -erf(-x)

它們與“Q 函數”密切相關,廣泛應用于通信工程學和正态分布內建。

// 累計分布(Cumulative distribution)

real normalDistributionFunction(real x) {

       return - 0.5 * erfc(-x * SQRT1_2);

}

real normalDistributionFunctionCompl(real x) {

       return -normalDistributionFunction(-x);

}

alias normalDistributionFunctionCompl Qfunction;

相反的累計正态分布函數,同時提供的還有被熟知的分位點函數(quantile function)。例如:如果考試分數依比例達60的平均分和12.5的标準方差,有學生分數會達到85或更多嗎?

import tango.io.Stdout;

import tango.math.ErrorFunction;

void main()

{

    // calculate how many standard deviations

    // we are from the mean.

    real x = (85-60)/12.5;

    Stdout.println(1-normalDistribution(x));

    // prints 0.0227501319481792, which tells us that 2.3% of students

    // will recieve a high distinction.

    // Now for the inverse calculation.

    // What mark is required to be in the top 1% of students?

    Stdout.println(60 + 12.5*normalDistributionInv(0.99));

    // prints: 89.0793484255105

math.GammaFunction

gamma, logGamma, sgnGamma

伽瑪函數是階乘函數包括分數和負數的一種歸納。在數理統計中這是極為重要的。 loggamma(x)與log(abs(gamma(x)))完全相等,但它有效地保留了更大的x的值。在C語言中,由于一些曆史原因把它們叫做tgamma和lgamma。sgnGamma(x)給出符号(x>0時gamma 通常是正的,但是符号有規律地變換為負x)。

beta, betaIncomplete, betaIncompleteInv

beta函數與伽瑪函數有密切關系。正如伽瑪函數gamma(x)類似x的階乘,beta(a,b)函數類似于二項式系數(binomial coefficient (a choose b)),但它對real型和整型(integers)數同樣有效。

不完全的beta函數betaIncomplete(a、b、x),是beta函數的一個部分,是最重要、廣泛應用于數理統計機率的計算。絕大多數基本的數理統計測試(如,學生的t,菲舍爾F測試( Fischer F tests))是通過不完全beta函數和它的逆運算(its inverse)來實作的。

gammaIncomplete, gammaIncompleteInv

不完全伽瑪函數(incomplete gamma function)在一些統計分布函數中用作低能級函數(low-level function ),其中最引人注目的是Chi-squared累積分布函數。

math.Bessel(貝塞耳函數)

貝塞爾函數是用來解決貝塞爾微分方程的:

x2 d2y/dx2 + x dy/dx + (x2 - a2)y = 0

這類方程最常見于圓柱或球形空間中的電磁波研究(studies of waves)(下一個管道,比如);是以,他們在電磁波的傳播學的計算中是很重要的,例如。有幾種類型的貝塞爾函數對應圓柱和球面坐标,不管x = 0時的值是finite ("first kind") 還是 infinite ("second kind")。傳統的數學中的這些函數的名字相當晦澀(J0 J1,Jn,Y0,Y1,Yn、)。在Tango中,他們被給予了輕松的更具描述性的名字[參考:WalterBrown]。

Value of a  First kind  Second kind

zero  cylBessel_J0  cylBessel_Y0

one  cylBessel_J1  cylBessel_Y1

integer n cylBessel_Jn  cylBessel_Yn

注:其它的貝塞爾( bessel)函數的變化不在目前的Tango實作中,但在将來可能添加到Tango的這個子產品裡。被修改的貝塞爾函數被重命名為cylBessel_I;第三種被修改的Bessel函數将被命名為cylBessel_K,球形的貝塞爾函數會被命名為sphBessel_j 和sphBessel_y

math.Elliptic(橢圓積分)

最初的橢圓積分出現在計算一個橢圓弧形長時,但通常适用于整個一類積分方程。任何此類方程可以轉換成這三種橢圓積分的結合。Tango的命名約定是基于數學應用軟體的。不完整的橢圓積分依賴振幅、φ和系數m 。完整的橢圓積分是φ= PI / 2時的特殊情形。

 Incomplete (不完整) Complete (完整)

First (第一) ellipticF  ellipticKComplete

Second (第二) ellipticE  ellipticEComplete

Third(第三)  ellipticPi*  ellipticPiComplete*

目前的實作受限于較低的精度。

在回顧

  對于準确性和快速代碼的一些一般性的建議:

  盡可能地利用零附近的相對高的精度。

  不要過分擔心除零錯誤。通常情況下,計算将依然是正确的。異常狀态可以用來檢測反常行為。

  正确使用NaN會非常有利于調試。

  你想要進行精确相等的測試使用feq忙天文超過= =,除非你想測試确切的平等。

  使用rndint和rndlong比用cast(int) 和 cast(long)好,除非平局決勝行為是很重要的。

  

附錄

 

 附錄1:Tango NaN結構

  “有效載荷”可以被儲存在NaN的尾部。一個位(最有意義的)必須區分NaN是安靜的還是在發信号。這就留下float的22位有效載荷,double的51位;80位real的62位;128位quad的111位。不幸的是,每當轉換到一個較小的類型的情況發生,位就會從有效載荷中丢失。在在Tango的NaN中, float的22位尾數是有效載荷整數的低22位。22..51位儲存在double的尾部,(80位real)位52 ..62儲存在real的尾部。

附錄2:數值算法目前還沒有包含在Tango中

值得注意的是,經典文本“數值分析方法庫”("Numerical Recipes",http://www.nrbook.com/a/),雖然是對數值算法進行廣泛的概述的優秀資源,但它某些部分包含一些錯誤的陳述,以及相關的源代碼有嚴重的著作權限制。

快速傅裡葉變換利用FFTW庫可以有效執行:http://www.fftw.org

 腳注:可移植性和80位實數的使用

為了可移植性,80位實數有時被認為應該犧牲。如果代碼是一個目标,輕便,認識到這一點很重要,盡管所有系統都能夠存儲資料和雙精度的浮動,并非所有的都能執行計算那些精度。特别的,x87浮點單元,現在幾乎所有的英特爾與AMD的機器的真實硬體支援80-bit浮點計算,在計算精度較低時包括一些仿真。奔騰III是英特爾的第一個原生支援雙精度浮點計算的處理器(需要使用SSE2指令集)。Tango庫盡可能使用用80位real型數,必要時優雅地退化到使用64位實數。

參考和延伸閱讀(References and Further Reading)

"What Every Computer Scientist Should Know About Floating-Point Arithmetic",

《每一個計算機科學家都應該知道什麼是浮點算法》

http://docs.sun.com/source/806-3568/ncg_goldberg.html

"A Proposal to Add Mathematical Special Functions to the C++ Standard Library", Walter E. Brown, 《建議增加數學特殊函數到c++标準庫》,沃爾特-布朗,

http://std.dkuug.dk/JTC1/SC22/WG21/docs/papers/2003/n1542.pdf. (Note that this proposal has been accepted into the C++0x TR1 standard,注意這個提案已經加入了c++ 0x TR1 标準).

"A Proposal to Add Mathematical Functions for Statistics to the C++ Standard Library", Paul A. Bristow,(《建議增加數理統計函數到c++标準庫》, Paul A. Bristow)

http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1069.pdf

"An Interview with the Old Man of Floating-Point: Reminiscences elicited from William Kahan by Charles Severance",

http://www.cs.berkeley.edu/~wkahan/ieee754status/754story.html

繼續閱讀