tango.math.GammaFunction

Implementation of the gamma and beta functions, and their integrals.

License:

BSD style: see license.txt

Authors:

Stephen L. Moshier (original C code). Conversion to D by Don Clugston
real MAXGAMMA [manifest]
The maximum value of x for which gamma(x) < real.infinity.
real sgnGamma(real x)
The sign of Γ(x).
Returns -1 if Γ(x) < 0, +1 if Γ(x) > 0, NAN if sign is indeterminate.
real gamma(real x)
The Gamma function, Γ(x)
Γ(x) is a generalisation of the factorial function to real and complex numbers. Like x!, Γ(x+1) = x*Γ(x).

Mathematically, if z.re > 0 then Γ(z) = 0 tz-1e-t dt

) ) ) ) ) ) )
Special Values
x Γ(x)
NAN NAN
±0.0 ±∞
integer > 0 (x-1)!
integer < 0 NAN
+∞ +∞
-∞ NAN
real logGamma(real x)
Natural logarithm of gamma function.
Returns the base e (2.718...) logarithm of the absolute value of the gamma function of the argument.

For reals, logGamma is equivalent to log(fabs(gamma(x))).

) ) ) )
Special Values
x logGamma(x)
NAN NAN
integer <= 0 +∞
±∞ +∞
real beta(real x, real y)
Beta function
The beta function is defined as

beta(x, y) = (Γ(x) Γ(y))/Γ(x + y)
real betaIncomplete(real aa, real bb, real xx)
Incomplete beta integral
Returns incomplete beta integral of the arguments, evaluated from zero to x. The regularized incomplete beta function is defined as

betaIncomplete(a, b, x) = Γ(a+b)/(Γ(a) Γ(b)) * 0x ta-1(1-t)b-1 dt

and is the same as the the cumulative distribution function.

The domain of definition is 0 <= x <= 1. In this implementation a and b are restricted to positive values. The integral from x to 1 may be obtained by the symmetry relation

betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x )

The integral is evaluated by a continued fraction expansion or, when b*x is small, by a power series.
real betaIncompleteInv(real aa, real bb, real yy0)
Inverse of incomplete beta integral
Given y, the function finds x such that

betaIncomplete(a, b, x) == y

Newton iterations or interval halving is used.
real gammaIncomplete(real a, real x)
real gammaIncompleteCompl(real a, real x)
Incomplete gamma integral and its complement
These functions are defined by

gammaIncomplete = ( 0x e-t ta-1 dt )/ Γ(a)

gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) = (x e-t ta-1 dt )/ Γ(a)

In this implementation both arguments must be positive. The integral is evaluated by either a power series or continued fraction expansion, depending on the relative values of a and x.
real gammaIncompleteComplInv(real a, real p)
Inverse of complemented incomplete gamma integral
Given a and y, the function finds x such that

gammaIncompleteCompl( a, x ) = p.

Starting with the approximate value x = a t3, where t = 1 - d - normalDistributionInv(p) sqrt(d), and d = 1/9a, the routine performs up to 10 Newton iterations to find the root of incompleteGammaCompl(a,x) - p = 0.
real digamma(real x)
Digamma function
The digamma function is the logarithmic derivative of the gamma function.

digamma(x) = d/dx logGamma(x)