Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

Log Gamma

概要

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

namespace boost{ namespace math{

template <class T>
calculated-result-type lgamma(T z);

template <class T, class Policy>
calculated-result-type lgamma(T z, const Policy&);

template <class T>
calculated-result-type lgamma(T z, int* sign);

template <class T, class Policy>
calculated-result-type lgamma(T z, int* sign, const Policy&);

}} // namespaces
說明

lgamma 函數 定義如下:

這個函數的第二種形式帶有一個指向一個整數的指針,如果這個指針不為空,那麼它就作為存儲tgamma(z)函數符號的輸出參數。

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

在庫的內部,這個函數有兩個高效的版本:一個完全通用的版本,運行速度慢一些但是更精確一些, 以及一個更高效的使用逼近方法的版本,當T的有效數字的個數與一個特定的 蘭克澤斯逼近對應時。實際上,任何一個你所遇到的內建的浮點類型都有一個適當的蘭克澤斯逼近 定義. 給定足夠的運算時間,也有可能使用libs/math/tools/lanczos_generator.cpp產生其它的蘭克澤斯逼近 .

函數的返回值的類型使用函數返回值推導法則 來確定:如果T是整型,那麼結果的類型是double,否則結果的類型是T.

精確性

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

注意:在lgamma函數正根附近的相對誤差很小,但是對於負值參數,lgamma函數有無限多個無理根: 在這個負根的附近僅可以保證較低的絕對誤差.

有效數字位數

平台和編譯器

階乘和半階乘

接近於0的值

接近於 1 或 2的值

接近於一個負極點的值

53

Win32 Visual C++ 8

峰值=0.88 均值=0.14

(GSL=33) (Cephes=1.5)

峰值=0.96 均值=0.46

(GSL=5.2) (Cephes=1.1)

峰值=0.86 均值=0.46

(GSL=1168) (Cephes~500000)

峰值=4.2 均值=1.3

(GSL=25) (Cephes=1.6)

64

Linux IA32 / GCC

峰值=1.9 均值=0.43

(GNU C Lib 峰值=1.7 均值=0.49)

峰值=1.4 均值=0.57

(GNU C Lib 峰值= 0.96 均值=0.54)

峰值=0.86 均值=0.35

(GNU C Lib 峰值=0.74 均值=0.26)

峰值=6.0 均值=1.8

(GNU C Lib 峰值=3.0 均值=0.86)

64

Linux IA64 / GCC

峰值=0.99 均值=0.12

(GNU C Lib 峰值 0)

峰值=1.2 均值=0.6

(GNU C Lib 峰值 0)

峰值=0.86 均值=0.16

(GNU C Lib 峰值 0)

峰值=2.3 均值=0.69

(GNU C Lib 峰值 0)

113

HPUX IA64, aCC A.06.06

峰值=0.96 均值=0.13

(HP-UX C Library 峰值 0)

峰值=0.99 均值=0.53

(HP-UX C Library 峰值 0)

峰值=0.9 均值=0.4

(HP-UX C Library 峰值 0)

峰值=3.0 均值=0.9

(HP-UX C Library 峰值 0)

測試

這些函數的主要測試包含與階乘的對數的比較,這些階乘可以獨立地計算到非常高的精度。

隨機測試同樣也使用了關鍵問題域(key problem areas)。

實現

通過合併針對於不完全γ函數的級數和連分數來實現這個函數的通用版本:

其中l 是一個任意的積分限制: 選擇l = max(10, a) 看起來效果很好. 對於負數z ,使用對數版本的反射方程:

對於精度已知的類型,使用蘭克澤斯逼近 ,一個 特性類(traits class ) boost::math::lanczos::lanczos_traits 用來將類型T與一個合適的逼近方法相關聯。對數版本的蘭克澤斯逼近 是:

其中 Le,g 是蘭克澤斯和, scaled by eg.

與前面一樣,當z < 0時使用這種方法。

當 z非常接近於 1 或 2, 對數版本的蘭克澤斯逼近 產生很大的消去誤差:實際上,對於非常接近於1或2的值,可能產生任意大的相對誤差(儘管絕對誤差很小).

對於精度高達113bit的類型(超過且包含128bit的long double 類型),在區間[1,2]和[2,3]上使用由JM 發明的根保留(root-preserving)有理逼近。超過區間[2,3]時使用的逼近形式如下:

lgamma(z) = (z-2)(z+1)(Y + R(z-2));

其中 Y 是一個常量, R(z-2)是有理逼近: 已經經過了優化,使得其絕對誤差與Y相比很小.些外,對於大於3的較小的z值可以使用下面的關係來進行參數簡化:

lgamma(z+1) = log(z) + lgamma(z);

在 區間[1,2]上,必須使用兩個逼近,一個是針對於較小的z值:

lgamma(z) = (z-1)(z-2)(Y + R(z-1));

Y是一個常量,R(z-1)針對於與Y相比的絕對誤差進行了優化。對於z>1.5 ,上面的形式不會收斂於一個極值解,但是下面的一個類似的形式可以做到這一點:

lgamma(z) = (2-z)(1-z)(Y + R(2-z));

最後,對於 z < 1 下面的遞推關係可以用來調整到 z > 1:

lgamma(z) = lgamma(z+1) - log(z);

注意:儘管這裡涉及到減法運算,但是並沒有產生消去誤差:當z從開始遞減,-log(z)項變為正值的速度比lgamma(z+1)變為負值的速度更快一些。因此,在這個特例中,有效數字是得到保留而不是消去。

對於那些有相對應的蘭克澤斯逼近 的類型,當前的解決方案如下:對於虛部,我們通過將冪項除以它在z=l處的值來平衡在 蘭克澤斯逼近 中的兩項, 然而我們將蘭克澤斯係數乘以相同的值。現在每一項都會在z = l處有1 這個值,且我們可以用過log1p的形式來重組這些冪項。類似地,如果我們從蘭克澤斯和中減去1 (在代數上,, 通過在z = 1將每一項減去這個值), 我們可以獲得一個用於傳遞給log1p的新值。關鍵地,所有的這些項都趨向於0 ,因為 z -> 1:

上面表達式中的Ck項與 蘭克澤斯逼近中的Ck項是一樣的。

可以在z = 2處進行類似的重排::


PrevUpHomeNext