Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

Digamma

概要

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

namespace boost{ namespace math{

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

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

}} // namespaces
說明

返回x的digamma 或 psi 函數. Digamma 定義為γ函數的對數導數。

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

這裡是這個函數的完全通用的版本:所有的這些版本都設置到一個指定的精度,最大精度為34位數字精度。

函數的返回值的類型使用函數返回值推導法則 來確定: 當T是整型的時候,返回值的類型是double,否則返回值的類型為T。

精確性

下面的表顯示了參數的不同定義域的峰值誤差(單位為10的-5次方)。注意,在那些由寬浮點類型來給定窄浮點類型的系統上會有有效的零誤差. 除非另外指定比此處的類型更窄的浮點類型,融計算結果具有有效的零誤差 .

有效數字位數

平台和編譯器

隨機正數值

接近於正根的值

接近於0的值

負值

53

Win32 Visual C++ 8

峰值=0.98 均值=0.36

峰值=0.99 均值=0.5

峰值=0.95 均值=0.5

峰值=214 均值=16

64

Linux IA32 / GCC

峰值=1.4 均值=0.4

峰值=1.3 均值=0.45

峰值=0.98 均值=0.35

峰值=180 均值=13

64

Linux IA64 / GCC

峰值=0.92 均值=0.4

峰值=1.3 均值=0.45

峰值=0.98 均值=0.4

峰值=180 均值=13

113

HPUX IA64, aCC A.06.06

峰值=0.9 均值=0.4

峰值=1.1 均值=0.5

峰值=0.99 均值=0.4

峰值=64 均值=6

就像上面顯示的那樣,對於正數參數,誤差率一般很低。對於負數參數有無限多個無理數根:相對誤差將會很大,而絕對誤差將會保持很低。

測試

混合的抽樣測試使用 functions.wolfram.com在線計算器產生的數據, 而隨機測試數據使用高精度的參考實現產生的數據(一個微分的 蘭克澤斯逼近,如下).

實現

實現被分為如下的定義域:

對於負值參數,使用下面的反射方程來使得x為正值:

digamma(1-x) = digamma(x) + pi/tan(pi*x);

對於在範圍[0,1]內的參數使用下面的迭代關係來將計算調整到範圍[1,2]內:

digamma(x) = digamma(x+1) - 1/x

對於在範圍 [1,2]內的參數,使用一個由 JM 發明的有理逼近 (如下).

對於在範圍 [2,BIG]內的參數使用下面的迭代關係來將計算調整到範圍[1,2]:

digamma(x+1) = digamma(x) + 1/x;

對於大於 BIG 的參數,可以使用下面的漸近展開:

然而,在一些項之後,這種展開有很多變化:到底有多少項依賴於x的大小。因此必須選定BIG的值,從而使得級數可以在一個在BIG值計算時足夠小以致於不會對計算結果產生任何影響的項上截斷。對於高達80-bit的實數將BIG的值選為10,對於高達128-bit的實數選擇BIG的值為20可以使得級數在一個合適的小數目的項後截斷且作為一個類似於在(1/x*x)的多項式內求值。

JM 發明的在區間[1 ,2]內的有理逼近導現如下.

首先,使用一個60項的 蘭克澤斯逼近來構造一個高精度的digamma逼近,形式如下:

其中P(x) 和 Q(x) 是源自於蘭克澤斯求和的無理數形式的多項式,而 P'(x) 和 Q'(x)是它們的一階導數。這個逼近的蘭克澤斯部分的理論逼近精度為~100個十進制數字。然而,上面求各式中的消去運算將會減少到99-(1/y),如果y是結果。這個逼近用來計算digamma的正根,並且發現與Cody所使用的25個數字一致(參見 Math. Comp. 27, 123-127 (1973) by Cody, Strecok 和 Thacher),也與Morris使用的35個數字一致(參見 TOMS Algorithm 708)。

類似地,一些抽樣測試與使用unctions.wolfram.com 超過40個數字的計算結果一致。這種足夠的精度可以確信逼近的精度可以達到double精度。達到128-bit long double精確度要求根的位置~70位數字,這種方法計算的數據是否滿足這個要求還不太清楚:困難在於獨立地獲取數據。

針對於絕對誤差的使用由 JM 發明的有理逼近來進行,形式如下:

digamma(x) = (x - X0)(Y + R(x - 1));

其中 X0 是digamma的正根, Y是一個常量, 而 R(x - 1) 是有理逼近. 注意,因為X0 是無理數,在X0中我們需要兩倍於x的數字位數來避免在減法中出現的消去誤差 (這裡假設x是一個精確值,如果不是這樣,那麼所有的假設都要取消 ). 這就意味著即使x是根捨入到最近的可表達的值,digamma(x)的結果也不會是0.


PrevUpHomeNext