std.mathspecial

pure nothrow @nogc @safe 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) = ∫₀^∞ t^z-1e^-t dt

Special Values

x	Γ(x)
NAN	NAN
±0.0	±∞
integer > 0	(x-1)!
integer < 0	NAN
+∞	+∞
-∞	NAN

pure nothrow @nogc @safe real logGamma(real x);

Natural logarithm of the gamma function, Γ(x)

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	+∞
±∞	+∞

pure nothrow @nogc @safe real sgnGamma(real x);

The sign of Γ(x).

Parameters:

real x

the argument of Γ

Returns:

-1 if Γ(x) < 0, +1 if Γ(x) > 0, and NAN if Γ(x) does not exist.

Note This function can be used in conjunction with logGamma to evaluate Γ(x) when gamma(x) is too large to be represented as a real.

Examples:

writeln(sgnGamma(10_000)); // 1

pure nothrow @nogc @safe real beta(real x, real y);

Beta function, B(x,y)

Mathematically, if x > 0 and y > 0 then B(x,y) = ∫₀¹t^x-1(l-t)^y-1dt. Through analytic continuation, it is extended to ℂ2 where it can be expressed in terms of Γ(z).

B(x,y) = Γ(x)Γ(y) / Γ(x+y).

This implementation restricts x and y to the set of real numbers.

Parameters:

real x	the first argument of B
real y	the second argument of B

Returns:

It returns B(x,y) if it can be computed, otherwise NAN.

Special Values

x	y	beta(x, y)
NAN	y	NAN
-∞	y	NAN
integer < 0	y	NAN
noninteger and x+y even ≤ 0	noninteger	-0
noninteger and x+y odd ≤ 0	noninteger	+0
+0	positive finite	+∞
+0	+∞	NAN
> 0	+∞	+0
-0	+0	NAN
-0	> 0	-∞
noninteger < 0, ⌈x⌉ odd	+∞	-∞
noninteger < 0, ⌈x⌉ even	+∞	+∞
noninteger < 0	±0	±∞

Since B(x,y) = B(y,x), if the table states that beta(x, y) is a special value, then beta(y, x) is one as well.

Examples:

writeln(beta(1, 2)); // 0.5

pure nothrow @nogc @safe real digamma(real x);

Digamma function, Ψ(x)

Ψ(x), is the logarithmic derivative of the gamma function, Γ(x).

Ψ(x) = d/_dx ln|Γ(x)| (the derivative of logGamma(x))

Parameters:

real x

the domain value

Returns:

It returns Ψ(x).

Special Values

x	digamma(x)
integer < 0	NAN
±0.0	∓∞
+∞	+∞
-∞	NAN
NAN	NAN

See Also:

logmdigamma, logmdigammaInverse.

Examples:

const euler = 0.57721_56649_01532_86060_65121L;

assert(isClose(digamma(1), -euler));
writeln(digamma(+0.)); // -real.infinity
writeln(digamma(-0.)); // +real.infinity
writeln(digamma(+real.infinity)); // +real.infinity
assert(isNaN(digamma(-1)));
assert(isNaN(digamma(-real.infinity)));

pure nothrow @nogc @safe real logmdigamma(real x);

Log Minus Digamma function

logmdigamma(x) = log(x) - digamma(x)

See Also:

digamma, logmdigammaInverse.

pure nothrow @nogc @safe real logmdigammaInverse(real x);

Inverse of the Log Minus Digamma function

Given y, the function finds x such log(x) - digamma(x) = y.

See Also:

logmdigamma, digamma.

pure nothrow @nogc @safe real betaIncomplete(real a, real b, real x);

Incomplete beta integral

Returns regularized 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) ) * ∫₀^x t^a-1(1-t)^b-1 dt

and is the same as the cumulative distribution function of the Beta distribution.

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.

pure nothrow @nogc @safe real betaIncompleteInverse(real a, real b, real y);

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.

pure nothrow @nogc @safe real gammaIncomplete(real a, real x);

pure nothrow @nogc @safe real gammaIncompleteCompl(real a, real x);

Incomplete gamma integral and its complement

These functions are defined by

gammaIncomplete = ( ∫₀^x e^-t t^a-1 dt )/ Γ(a)

gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) = (∫_x^∞ e^-t t^a-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.

pure nothrow @nogc @safe real gammaIncompleteComplInverse(real a, real p);

Inverse of complemented incomplete gamma integral

Given a and p, the function finds x such that

gammaIncompleteCompl( a, x ) = p.

pure nothrow @nogc @safe real erf(real x);

Error function

The integral is

erf(x) = 2/ √(π) ∫₀^x exp( - t²) dt

The magnitude of x is limited to about 106.56 for IEEE 80-bit arithmetic; 1 or -1 is returned outside this range.

pure nothrow @nogc @safe real erfc(real x);

Complementary error function

erfc(x) = 1 - erf(x) = 2/ √(π) ∫_x^∞ exp( - t²) dt

This function has high relative accuracy for values of x far from zero. (For values near zero, use erf(x)).

pure nothrow @nogc @safe real normalDistribution(real x);

Standard normal distribution function.

The normal (or Gaussian, or bell-shaped) distribution is defined as:

normalDist(x) = 1/√(2π) ∫_-∞^x exp( - t²/2) dt = 0.5 + 0.5 * erf(x/sqrt(2)) = 0.5 * erfc(- x/sqrt(2))

To maintain accuracy at values of x near 1.0, use normalDistribution(x) = 1.0 - normalDistribution(-x).

References http://www.netlib.org/cephes/ldoubdoc.html, G. Marsaglia, "Evaluating the Normal Distribution", Journal of Statistical Software 11, (July 2004).

pure nothrow @nogc @safe real normalDistributionInverse(real p);

Inverse of Standard normal distribution function

Returns the argument, x, for which the area under the Normal probability density function (integrated from minus infinity to x) is equal to p.

Note This function is only implemented to 80 bit precision.