動機(Motivation)

為什麼讓γ及類似γ的函數基於蘭克澤斯逼近？

首先，我必須說明的是在實數域上的函數使用蘭克澤斯逼近(參考Wikipedia 或 Mathworld) ，對於其它傳統的方法，比如Stirling's approximation，並沒有提供很明顯的優點。Pugh 為多種可用的方法進行了一個擴展的比較，並發現在誤差和複雜性方面它們都是很類似的。然而，蘭克澤斯逼近卻有一些特性使得它值得進一步的考慮：

這兩個屬性的組合使得這個逼近方法非常具有吸引力：斯特林( Stirling approximation)逼近對於較大的z值很精確，並且與蘭克澤斯逼近有一些相同的分析性質，但是不能在整個z的範圍上使用。

作為一個簡單的例子，考慮兩個γ函數的比：可以使用 lgamma 函數來計算：

exp(lgamma(a) - lgamma(b));

然而，即使 lgamma 是一致地精確到0.5ulp，上面的方法的最嚴重的相對誤差，可以寫作：

Erel > a * log(a)/2 + b * log(b)/2

對於小的a 和 b ，這並不是一個問題，但以另外一種方式來看待這種關係：每次以10為因子來增大 a 和 b 的值，那麼至少有一個十進制數字的精度丟失。

相反，通過組合相似的冪項為一個蘭克澤斯逼近，對於較小的a 和 b，這些誤差可以消除，並且對於非常大的a和b，可以使誤差在可控制的範圍之內。當然，計算次數巨大 的冪本來就是一個很艱難的問題，但即使是這樣，蘭克澤斯分析組合可以在一個合理的結果和一個無用值之間產生差別。參考β 函數的實現說明來瞭解實際中使用這個方法的例子，不完全gamma_p gamma 和 β 函數使用類似的冪項的分析組合，將被大的冪項除的γ和β函數組合成一個單獨的(更簡單的)表達式。

逼近

對數的蘭克澤斯逼近由下面的方程給定：

其中 S_g(z) 是一個無限求和，對於 所有的 z > 0收斂，且g 是一個控制這個求和式中的項的「形狀」的任何的參數：

單獨的係數由閉型(closed form)定義為：

然而，在計算上長階乘和下降階乘的過程中，以這種形式的求和方式會導致數值不穩定性(我們可以有效地乘以一系列的非常接近於1的值，所以捨入誤差可以很快進行統計。)。

蘭克澤斯逼近通常寫作部分階乘的形式，前導常量被求和式中的係數所吸收：

其中：

參數g 是一個任意選定的常量，且N 是在「蘭克澤斯和」中任意選定的計算的項的數目。

	注意
	一些作者選擇定義從1到N定義和式，因此有N＋1個係數。這會使下面的討論和代碼發生混淆(因為C++處理關開的數級範圍而不是求和式的閉全範圍)。這種約定形式與Godfrey一致，但是與Pugh不一致，因此在涉及到這些術語的時候要注意。

計算係數(Computing the Coefficients)

係數 C0..CN-1 在精度較高時從N 和g 中進行計算，然後存儲為程序的一部分。計算係數使用Godfrey方法;修定常量被存儲在一個列向量P中，那麼：

P = B D C F

其中 B是一個 NxN 矩陣：

D 是一個 NxN 矩陣：

C 是一個 NxN 矩陣：

F是一個有N個元素的列向量：

注意矩陣 B, D 和 C 包含所有的項數並只依賴於N，這種結果必須首先計算然而乘以F作為最後一步。

選擇正確的參數(Choosing the Right Parameters)

這裡的技巧是選擇N 和 g 來給定預期的精度：為g選擇一個較小的值將產生一個嚴格收斂的級數，但收斂很慢。為g選擇一個較大的值，導致級數的項數非常大，並且對於開始的g-1 項是多種多樣的，然後突然收斂。

Pugh 對於在範圍 1 <= N <= 60中的N值推導出了優化的g值：不幸的是，在實際中選擇這些值會導致在蘭克澤斯求和中發生消去錯誤，錯誤中的最大項大約比結果大1000倍。除非這些計算可以使用一些guard數字並將這些係數以高於結果精度的精度存儲，否則這些優化的值在實際中是沒有用的。這些值最好為使用double精度的算術來計算float精度。

由Godfrey 描述的另一種方法是：對於一個給定p個數字的浮點類型，在N和p的參數空間上進行足夠多的查找，來找出最佳的N和g的組合。重複這種操作，對於double精度的算術找到最好的逼近 (接近於Godfrey 找到的一個)，但是對於80-bit或128-bit的long double類型不能找到真正好的逼近。此外，獲得的逼近趨向於對用於測試這些實現的較小的z值更優化一些。用巨大的參數來計算函數的比被發現會因為對蘭克澤斯級數的截斷而產生誤差。

Pugh 鑒定出了所有的逼近的理論誤差的最小值的位置，但不幸的是，只發佈了這些極小值中的最大的那一部分。然而，他觀測到：極值非常接近蘭克澤斯逼近中的第一個被忽略的項。這些位置很容易定位，儘管需要大量的計算時間。這些"sweet spots" 只需要計算一次，製表，然而對於某些固定精度的類型進行所需要的精度的逼近時進行查找。

不幸的是，沿著這種方式不能為128-bit的long double類型找到真正好的逼近，並且為64-bit和80-bit的實數類型查找逼近需要大量的係數項。這有兩個競爭的因素：更高的精度要求巨大的g值，但是避免消去錯誤要求一個小的g值。

在這裡，注意蘭克澤斯逼近可以轉換為有理形式(兩個多項式的比，使用多項式算術從部分階乘的形式中得到)，並且這樣做會改變係數，使得它們都為正數。這意味著在有理形式中的和可以在沒有消去錯誤的情況下計算得到，即使對於一個給定的N將係數的個數加倍。重複查找"sweet spots"，這一次使用有理函數的形式來計算蘭克澤斯逼近，並且只測試這些理論誤差與正大測試的類型的機器精度更小的"sweet spots"，對於所有的測試的類型都產生了很好的逼近。這些最佳值與Pugh 報告的最佳情況非常吻合(對於一個給定的精度，N 比 Pugh 的報告值要稍大一些，g 比Pugh 報告的值要稍小一些)，儘管轉換為有理函數形式將係數個數增加了兩倍，但應當指出的是它們中有一半是整數(因此只需要更少的存儲空間)並且將會要求使用較小的N值，所以只需要使用更少的浮點運算操作。

當以一個固定精度來計算時，下面的表中顯示了N 和 g 的最佳值。這些值應當在配合運算時使用：對於106-bit的機器沒有最佳值 (Darwin long doubles & NTL quad_float)，並且對g值的近一步優化也是可能的。這個表中的誤差是將一個無限蘭克澤斯逼近截斷到N項時的估計誤差。誤差從最開始的5個被忽略的項計算得來－並已知是非常悲觀的估計－當估計的截斷誤差與問題中所討論的類型的機器精度幾乎完全相符時，N和g的最佳組合就出現了。

最後請注意，通過將一個因子從exp(g)的分母中刪除可以將蘭克澤斯逼近寫為下面的形式，然後用exp(g)除以所有的係數：

這種形式對於計算 lgamma函數更方便，但是對於 gamma 函數，除以e值會導致將一個可能完全精確的值變為一個不精確的值：這會在輸入為精確值的情況下降低精度，因此這種形式沒有用在 gamma 函數中。

參考資料