Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

不完全γ函數

概要

#include <boost/math/special_functions/gamma.hpp>

namespace boost{ namespace math{

template <class T1, class T2>
calculated-result-type gamma_p(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type gamma_p(T1 a, T2 z, const Policy&);

template <class T1, class T2>
calculated-result-type gamma_q(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type gamma_q(T1 a, T2 z, const Policy&);

template <class T1, class T2>
calculated-result-type tgamma_lower(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type tgamma_lower(T1 a, T2 z, const Policy&);

template <class T1, class T2>
calculated-result-type tgamma(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type tgamma(T1 a, T2 z, const Policy&);

}} // namespaces
說明

有四個不完全γ函數:兩個是規範形式(也被稱作規範不完全函數),返回範圍在[0, 1]區間內的值,而另外兩個是不規則形式,返回值的範圍為 [0, Γ(a)]。對於統計應用感興趣的用戶應當使用 規範形式 (gamma_p 和 gamma_q).

所有的這些函數要求a > 0z >= 0, 否則返回定義域錯誤.

最後一個策略 參數是可選的並且可以用來控制函數的行為: 如何處理錯誤, 使用哪種層次的精度等等. 參見策略文檔瞭解更多信息。

當T1和T2是不同類型的時候,函數返回值的類型使用返回值推導法則:來確定, 否則返回值的類型是T1.

template <class T1, class T2>
calculated-result-type gamma_p(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type gamma_p(T1 a, T2 z, const Policy&);

返回a和z的規範的下半部不完全γ函數:

圍繞在點z==a處,函數很快地從0變化到1:

template <class T1, class T2>
calculated-result-type gamma_q(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type gamma_q(T1 a, T2 z, const Policy&);

返回a和z的規範的上半部不完全γ函數:

圍繞在點z==a處,函數很快地從1變化到0:

template <class T1, class T2>
calculated-result-type tgamma_lower(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type tgamma_lower(T1 a, T2 z, const Policy&);

返回完全的(非規範化) a和z的下半部不完全γ函數:

template <class T1, class T2>
calculated-result-type tgamma(T1 a, T2 z);

template <class T1, class T2, class Policy>
calculated-result-type tgamma(T1 a, T2 z, const Policy&);

返回完全的(非規範化) a和z的上半部不完全γ函數:

精確度

下面的圖像顯示了在a和z的不同的定義域上的函數的相對誤差的峰值和均值,以及與GSL-1.9 庫和Cephes 庫的比較.注意:在一些平台上使用窄浮點類型來表示寬浮點類型的唯一結果具有高效零誤差.

注意:當 a增大時,誤差也在增大 。

同樣也要注意:對於80bit 精度和120-bit精度的高誤差具有一些迷惑性:在64-bit期待為0的結果在更大的指數範圍long double 可能不為零,但是會非常小,這些函數也因此反映了為這些類型構造的測試的更極限的屬性。

所有的數據都是以epsilon(10的-5次方)為單位。

表8.函數 gamma_p(a,z)的誤差

有效數字位數

平台和編譯器

0.5 < a < 100

0.01*a < z < 100*a

1x10-12 < a < 5x10-2

0.01*a < z < 100*a

1e-6 < a < 1.7x106

1 < z < 100*a

53

Win32, Visual C++ 8

峰值=36 均值=9.1

(GSL 峰值=342 均值=46)

(Cephes 峰值=491 均值=102)

峰值=4.5 均值=1.4

(GSL 峰值=4.8 均值=0.76)

(Cephes 峰值=21 均值=5.6)

峰值=244 均值=21

(GSL 峰值=1022 均值=1054)

(Cephes 峰值~8x106 均值~7x104)

64

RedHat Linux IA32, gcc-3.3

峰值=241 均值=36

峰值=4.7 均值=1.5

峰值~30,220 均值=1929

64

Redhat Linux IA64, gcc-3.4

峰值=41 均值=10

峰值=4.7 均值=1.4

峰值~30,790 均值=1864

113

HPUX IA64, aCC A.06.06

峰值=40.2 均值=10.2

峰值=5 均值=1.6

峰值=5,476 均值=440


表9.函數 gamma_q(a,z)的誤差

有效數字位數

平台和編譯器

0.5 < a < 100

0.01*a < z < 100*a

1x10-12 < a < 5x10-2

0.01*a < z < 100*a

1x10-6 < a < 1.7x106

1 < z < 100*a

53

Win32, Visual C++ 8

峰值=28.3 均值=7.2

(GSL 峰值=201 均值=13)

(Cephes 峰值=556 均值=97)

峰值=4.8 均值=1.6

(GSL 峰值~1.3x1010 均值=1x10+9)

(Cephes 峰值~3x1011 均值=4x1010)

峰值=469 均值=33

(GSL 峰值=27,050 均值=2159)

(Cephes 峰值~8x106 均值~7x105)

64

RedHat Linux IA32, gcc-3.3

峰值=280 均值=33

峰值=4.1 均值=1.6

峰值=11,490 均值=732

64

Redhat Linux IA64, gcc-3.4

峰值=32 均值=9.4

峰值=4.7 均值=1.5

峰值=6815 均值=414

113

HPUX IA64, aCC A.06.06

峰值=37 均值=10

峰值=11.2 均值=2.0

峰值=4,999 均值=298


Table?0.Function tgamma_lower(a,z)誤差

有效數字位數

平台和編譯器

0.5 < a < 100

0.01*a < z < 100*a

1x10-12 < a < 5x10-2

0.01*a < z < 100*a

53

Win32, Visual C++ 8

峰值=5.5 均值=1.4

峰值=3.6 均值=0.78

64

RedHat Linux IA32, gcc-3.3

峰值=402 均值=79

峰值=3.4 均值=0.8

64

Redhat Linux IA64, gcc-3.4

峰值=6.8 均值=1.4

峰值=3.4 均值=0.78

113

HPUX IA64, aCC A.06.06

峰值=6.1 均值=1.8

峰值=3.7 均值=0.89


表10.函數 tgamma(a,z)的誤差

有效數字位數

平台和編譯器

0.5 < a < 100

and

0.01*a < z < 100*a

1x10-12 < a < 5x10-2

and

0.01*a < z < 100*a

53

Win32, Visual C++ 8

峰值=5.9 均值=1.5

峰值=1.8 均值=0.6

64

RedHat Linux IA32, gcc-3.3

峰值=596 均值=116

峰值=3.2 均值=0.84

64

Redhat Linux IA64, gcc-3.4.4

峰值=40.2 均值=2.5

峰值=3.2 均值=0.8

113

HPUX IA64, aCC A.06.06

峰值=364 均值=17.6

峰值=12.7 均值=1.8


測試

有兩種類型的測試:抽樣測試的比較數據基於這個庫的實現並使用數學 世界在線計算器 用來進行一些基本的"合理性檢查(sanity check)".精度的測試數據使用非常高精度的測試數據 (使用精度為1000-bit的 NTL's RR 類) ,使用一個60項的蘭克澤斯逼近產生非常高精度的數據,同時有一些但不是所有的特殊情況都被禁止了。這還不太令人滿意:應當使用一個獨立的方法,但是缺少這一些獨立的方法。我們甚至不能故意單純使用一個不對特殊情況進行處理的庫實現,因為對於較小的a 和 z ,勒讓德連分數是不穩定的。

實現

這四個函數共同使用一個通用的實現,因為它們通過下面的關係聯繫起來:

1)

2)

3)

不完全γ函數的下半部使用它的級數表示來計算:

4)

或者當 x > a 且 x > 1.1時從上半部積分中減去Γ(a) 或 1.

上半部積分使用勒讓德連分數表示來計算:

5)

x < a 時從下半部積分中減去Γ(a) 或 1.

對於x < 1.1,計算上半部積分有一些複雜,因為在這個區域內的連分數表示並不是連續的。然而還有另外一個級數表示可用於下半部積分:

6)

這就可以通過重組來計算上半部積分:

7)

參看powm1tgamma1pm1瞭解這些函數的實現.注意: tgamma1pm1 的精度交叉於35個數字或與與T類型相關的蘭克澤斯逼近 相關- 如果有的話 - 無論哪一個更大一些. 這就為這個區域內的精度強加了一個類似的限制。

對於x < 1.1 結果~0.5的交叉點不會在x~y時出現。使用x * 1.1 < a 作為當 0.5 < x <= 1.1 時的交叉準則保持最大的計算值在0.6附近。類似的,對於x <= 0.5 使用-0.4 / log(x) < a作為交叉準則保持最大的計算值在0.7附近 (不管是上半部積分還是下半部積分).

當a是一個整數或半整數的時候使用了兩個特殊情況,且上面列舉的交叉情況表明我們應當計算上半部積分Q。如果a是一個在範圍1 <= a < 30 中的值,那麼使用下面的有限求和:

9)

而對於在範圍0.5 <= a < 30 內的值,使用下面的有限求和:

10)

它們比連分數法更穩定也更高效一些。

當參數a很大時,且x~a,那麼級數(4)和連分數(5)收斂得很慢,在這個區域使用坦普爾展開:

11)

12)

13)

14)

被截斷為固定數目的項-用於產生一個指定的精確度-且作為一個多項式的式項式來求值。有高達128-bit精確度的函數版本:對於那些要求更高精度的類型不使用上面的展開方法。係數Ckn 由Temme給定的遞推關係預先計算,使用這些展開方法為:

(a > 20) && (a < 200) && fabs(x-a)/a < 0.4

以及:

(a > 200) && (fabs(x-a)/a < 4.5/sqrt(a))

後面的一個範圍對於所有的高達128-bit精度的long double的類型都是合法的,並且被設計用來確保計算結果大於10-6,第一個範圍只用於高達80-bit的long double類型。這些定義域比Temme 或 Didonato 和 Morris推薦的定義域都要窄。然而,使用一個範圍導致更大的且更不精確(也就是說,計算的結果)的值被傳遞給函數exp和函數erfc,從而導致更大的誤差率。換句話說,在效率和誤差之間有一個權衡。當前的當前為保持被(4)和(5)要求的項數在double精度時不超過~20。

對於規範的不完全函數,計算前導冪係數是這個函數的精度的中心。對於小的a和x合併冪集項使用蘭克澤斯逼近 給出最大的精度:

15)

結果這會產生下溢和上溢,然後這個指數可以消去一個因子 la 並放在冪項裡面

當a和x很大時,我們會得到一個非常大的指數,其基數接近於1:由冪函數來計算不會產生精確的結果,而簡單地取對數導致消去誤差。最嚴重的誤差可以使用下面的方法來避免:

16)

a-x 很小而 a 和 x 很大 .這裡仍然會有一個減法操作因此,會產生一些消去誤差 - 因為項很小,所以絕對誤差將會很小,是絕對誤差而不是相對誤差在函數exp中讀數。注意,對於足夠大的a和x,你仍然會得到誤差,雖然與其它方法相比這會延遲這些誤差不避免的誤差。在這裡使用log(1+x)-x 是受到Temme 的啟發(參看下面的參考資料).

參考資料

PrevUpHomeNext