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);

Regularized incomplete beta function I_x(a,b)

Mathematically, if a and b are positive real numbers, and 0 ≤ x ≤ 1, then I_x(a,b) = ∫₀^xt^a-1(1-t)^b-1dt/B(a,b) where B is the beta function. It is also the cumulative distribution function of the beta distribution.

betaIncomplete(a, b, x) evaluates I_x(a,b).

Parameters:

real a	the first argument of B, must be positive
real b	the second argument of B, must be positive
real x	the fraction of integration completion from below, 0 ≤ x ≤ 1

Returns:

It returns I_x(a,b), an element of [0,1].

Special Values

a	b	x	betaIncomplete(a, b, x)
negative	b	x	NAN
a	negative	x	NAN
a	b	< 0	NAN
a	b	> 1	NAN
+0	+0	(0,1)	NAN
∞	∞	(0,1)	NAN

If one or more of the input parameters are NAN, the one with the largest payload is returned. For equal payloads but with possibly different signs, the order of preference is x, a, b.

Note The integral is evaluated by a continued fraction expansion or, when b * x is small, by a power series.

See Also:

beta betaIncompleteCompl

Examples:

writeln(betaIncomplete(1, 1, .5)); // .5
writeln(betaIncomplete(+0., +0., 0)); // 0
assert(isNaN(betaIncomplete(+0., +0., .5)));
assert(isNaN(betaIncomplete(real.infinity, real.infinity, .5)));
writeln(betaIncomplete(real.infinity, real.infinity, 1)); // 1
assert(betaIncomplete(NaN(0x1), 1, NaN(0x2)) is NaN(0x2));
assert(betaIncomplete(1, NaN(0x3), -NaN(0x3)) is -NaN(0x3));

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

Regularized incomplete beta function complement IC_x(a,b)

Mathematically, if a > 0, b > 0, and 0 ≤ x ≤ 1, then IC_x(a,b) = ∫_x¹t^a-1(1-t)^b-1dt/B(a,b) where B is the beta function. It is also the complement of the cumulative distribution function of the beta distribution. It can be shown that IC_x(a,b) = I_1-x(b,a).

betaIncompleteCompl(a, b, x) evaluates IC_x(a,b).

Parameters:

real a	the first argument of B, must be positive
real b	the second argument of B, must be positive
real x	the fraction of integration completion from above, 0 ≤ x ≤ 1

Returns:

It returns IC_x(a,b), an element of [0,1].

Special Values

a	b	x	betaIncompleteCompl(a, b, x)
negative	b	x	NAN
a	negative	x	NAN
a	b	< 0	NAN
a	b	> 1	NAN
+0	+0	(0,1)	NAN
∞	∞	(0,1)	NAN

If one or more of the input parameters are NAN, the one with the largest payload is returned. For equal payloads but with possibly different signs, the order of preference is x, a, b.

See Also:

beta betaIncomplete

Examples:

writeln(betaIncompleteCompl(.1, .2, 0)); // betaIncomplete(.2, .1, 1)

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);

Regularized lower incomplete gamma function P(a,x)

Mathematically, P(a,x) = γ(a,x)/Γ(a), where γ(a,x) is the lower incomplete gamma function, γ(a,x) = ∫₀^xt^a-1e^-tdt, a > 0, and x ≥ 0.

Parameters:

real a	the shape parameter, must be positive
real x	the fraction of integration completion from below, must be non-negative

Returns:

It returns P(a,x), an element of [0,1].

Special Values

a	x	gammaIncomplete(a, x)
negative		NAN
	< 0	NAN
positive	0	0
positive	∞	1
+0	> 0	1
∞	(0, ∞)	0

See Also:

gamma and gammaIncompleteCompl

Examples:

assert(isClose(gammaIncomplete(1, 1), 1 - 1/E));
writeln(gammaIncomplete(1, 0)); // 0
writeln(gammaIncomplete(1, real.infinity)); // 1
writeln(gammaIncomplete(+0., 1)); // 1
writeln(gammaIncomplete(real.infinity, 1)); // 0

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

Regularized upper incomplete gamma function Q(a,x)

Mathematically, Q(a,x) = Γ(a,x)/Γ(a), where Γ(a,x) is the upper incomplete gamma function, Γ(a,x) = ∫_x^∞t^a-1e^-tdt, a > 0, and x ≥ 0. Note that P(a,x) + Q(a,x) = 1 or Q(a,x) = 1 - P(a,x), so Q is the complement of P.

Parameters:

real a	the shape parameter, must be positive
real x	the fraction of integration completion from above, must be non-negative

Returns:

It returns Q(a,x), an element of [0,1].

Special Values

a	x	gammaIncompleteCompl(a, x)
negative		NAN
	< 0	NAN
positive	0	1
positive	∞	0
+0	> 0	0
∞	(0, ∞)	1

See Also:

gamma and gammaIncomplete

Examples:

assert(isClose(gammaIncompleteCompl(2, 1), 2/E));
writeln(gammaIncompleteCompl(1, 0)); // 1
writeln(gammaIncompleteCompl(1, real.infinity)); // 0
writeln(gammaIncompleteCompl(+0., 1)); // 0
writeln(gammaIncompleteCompl(real.infinity, 1)); // 1
assert(isClose(gammaIncompleteCompl(1, 2), 1-gammaIncomplete(1, 2)));

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

Inverse regularized upper incomplete gamma function Q-1(a,p) with respect to p

Given a and p, the function finds x such that p = Q(a,x).

Parameters:

real a	the shape parameter, must be positive
real p	Q(a,x), must be in the interval [0,1]

Returns:

It returns x, a value ≥ 0

Special Values

a	p	gammaIncompleteComplInverse(a, p)
negative		NAN
	< 0	NAN
	> 1	NAN
+0	< 1	NAN
∞	> 0	NAN
> 0	0	∞
< ∞	1	0

See Also:

gammaIncompleteCompl

Examples:

const a = 2, p = 0.5L;
assert(isClose(gammaIncompleteComplInverse(a, gammaIncompleteCompl(a, p)), p));

assert(isClose(gammaIncompleteComplInverse(1, 1/E), 1));
assert(isNaN(gammaIncompleteComplInverse(+0.0L, 0.1)));
assert(isNaN(gammaIncompleteComplInverse(real.infinity, 0.2)));
assert(gammaIncompleteComplInverse(3, 0) is real.infinity);
writeln(gammaIncompleteComplInverse(4, 1)); // 0

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.