1 // Written in the D programming language.
4 * Mathematical Special Functions
6 * The technical term 'Special Functions' includes several families of
7 * transcendental functions, which have important applications in particular
8 * branches of mathematics and physics.
10 * The gamma and related functions, and the error function are crucial for
11 * mathematical statistics.
12 * The Bessel and related functions arise in problems involving wave propagation
13 * (especially in optics).
14 * Other major categories of special functions include the elliptic integrals
15 * (related to the arc length of an ellipse), and the hypergeometric functions.
18 * Many more functions will be added to this module.
19 * The naming convention for the distribution functions (gammaIncomplete, etc)
20 * is not yet finalized and will probably change.
23 * TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
24 * <caption>Special Values</caption>
26 * SVH = $(TR $(TH $1) $(TH $2))
27 * SV = $(TR $(TD $1) $(TD $2))
30 * SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
34 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>)
35 * POWER = $1<sup>$2</sup>
36 * SUB = $1<sub>$2</sub>
37 * BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>)
38 * CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG ))
41 * PLUSMNINF = ±∞
49 * Copyright: Based on the CEPHES math library, which is
50 * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
51 * License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
52 * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston
53 * Source: $(PHOBOSSRC std/_mathspecial.d)
55 module std
.mathspecial
;
56 import std
.internal
.math
.errorfunction
;
57 import std
.internal
.math
.gammafunction
;
58 public import std
.math
;
60 /* ***********************************************
61 * GAMMA AND RELATED FUNCTIONS *
62 * ***********************************************/
69 /** The Gamma function, $(GAMMA)(x)
71 * $(GAMMA)(x) is a generalisation of the factorial function
72 * to real and complex numbers.
73 * Like x!, $(GAMMA)(x+1) = x * $(GAMMA)(x).
75 * Mathematically, if z.re > 0 then
76 * $(GAMMA)(z) = $(INTEGRATE 0, $(INFIN)) $(POWER t, z-1)$(POWER e, -t) dt
79 * $(SVH x, $(GAMMA)(x) )
80 * $(SV $(NAN), $(NAN) )
81 * $(SV $(PLUSMN)0.0, $(PLUSMNINF))
82 * $(SV integer > 0, (x-1)! )
83 * $(SV integer < 0, $(NAN) )
84 * $(SV +$(INFIN), +$(INFIN) )
85 * $(SV -$(INFIN), $(NAN) )
90 return std
.internal
.math
.gammafunction
.gamma(x
);
93 /** Natural logarithm of the gamma function, $(GAMMA)(x)
95 * Returns the base e (2.718...) logarithm of the absolute
96 * value of the gamma function of the argument.
98 * For reals, logGamma is equivalent to log(fabs(gamma(x))).
101 * $(SVH x, logGamma(x) )
102 * $(SV $(NAN), $(NAN) )
103 * $(SV integer <= 0, +$(INFIN) )
104 * $(SV $(PLUSMNINF), +$(INFIN) )
107 real logGamma(real x
)
109 return std
.internal
.math
.gammafunction
.logGamma(x
);
112 /** The sign of $(GAMMA)(x).
114 * Returns -1 if $(GAMMA)(x) < 0, +1 if $(GAMMA)(x) > 0,
115 * $(NAN) if sign is indeterminate.
117 * Note that this function can be used in conjunction with logGamma(x) to
118 * evaluate gamma for very large values of x.
120 real sgnGamma(real x
)
122 /* Author: Don Clugston. */
123 if (isNaN(x
)) return x
;
124 if (x
> 0) return 1.0;
125 if (x
< -1/real.epsilon
)
127 // Large negatives lose all precision
133 return x
== 0 ?
copysign(1, x
) : real.nan
;
135 return n
& 1 ?
1.0 : -1.0;
140 assert(sgnGamma(5.0) == 1.0);
141 assert(isNaN(sgnGamma(-3.0)));
142 assert(sgnGamma(-0.1) == -1.0);
143 assert(sgnGamma(-55.1) == 1.0);
144 assert(isNaN(sgnGamma(-real.infinity
)));
145 assert(isIdentical(sgnGamma(NaN(0xABC)), NaN(0xABC)));
150 * The beta function is defined as
152 * beta(x, y) = ($(GAMMA)(x) * $(GAMMA)(y)) / $(GAMMA)(x + y)
154 real beta(real x
, real y
)
158 return exp(logGamma(x
) + logGamma(y
) - logGamma(x
+y
));
159 } else return gamma(x
) * gamma(y
) / gamma(x
+y
);
164 assert(isIdentical(beta(NaN(0xABC), 4), NaN(0xABC)));
165 assert(isIdentical(beta(2, NaN(0xABC)), NaN(0xABC)));
170 * The digamma function is the logarithmic derivative of the gamma function.
172 * digamma(x) = d/dx logGamma(x)
174 * See_Also: $(LREF logmdigamma), $(LREF logmdigammaInverse).
178 return std
.internal
.math
.gammafunction
.digamma(x
);
181 /** Log Minus Digamma function
183 * logmdigamma(x) = log(x) - digamma(x)
185 * See_Also: $(LREF digamma), $(LREF logmdigammaInverse).
187 real logmdigamma(real x
)
189 return std
.internal
.math
.gammafunction
.logmdigamma(x
);
192 /** Inverse of the Log Minus Digamma function
194 * Given y, the function finds x such log(x) - digamma(x) = y.
196 * See_Also: $(LREF logmdigamma), $(LREF digamma).
198 real logmdigammaInverse(real x
)
200 return std
.internal
.math
.gammafunction
.logmdigammaInverse(x
);
203 /** Incomplete beta integral
205 * Returns incomplete beta integral of the arguments, evaluated
206 * from zero to x. The regularized incomplete beta function is defined as
208 * betaIncomplete(a, b, x) = $(GAMMA)(a + b) / ( $(GAMMA)(a) $(GAMMA)(b) ) *
209 * $(INTEGRATE 0, x) $(POWER t, a-1)$(POWER (1-t), b-1) dt
211 * and is the same as the the cumulative distribution function.
213 * The domain of definition is 0 <= x <= 1. In this
214 * implementation a and b are restricted to positive values.
215 * The integral from x to 1 may be obtained by the symmetry
218 * betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x )
220 * The integral is evaluated by a continued fraction expansion
221 * or, when b * x is small, by a power series.
223 real betaIncomplete(real a
, real b
, real x
)
225 return std
.internal
.math
.gammafunction
.betaIncomplete(a
, b
, x
);
228 /** Inverse of incomplete beta integral
230 * Given y, the function finds x such that
232 * betaIncomplete(a, b, x) == y
234 * Newton iterations or interval halving is used.
236 real betaIncompleteInverse(real a
, real b
, real y
)
238 return std
.internal
.math
.gammafunction
.betaIncompleteInv(a
, b
, y
);
241 /** Incomplete gamma integral and its complement
243 * These functions are defined by
245 * gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
247 * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x)
248 * = ($(INTEGRATE x, $(INFIN)) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
250 * In this implementation both arguments must be positive.
251 * The integral is evaluated by either a power series or
252 * continued fraction expansion, depending on the relative
255 real gammaIncomplete(real a
, real x
)
261 return std
.internal
.math
.gammafunction
.gammaIncomplete(a
, x
);
265 real gammaIncompleteCompl(real a
, real x
)
271 return std
.internal
.math
.gammafunction
.gammaIncompleteCompl(a
, x
);
274 /** Inverse of complemented incomplete gamma integral
276 * Given a and p, the function finds x such that
278 * gammaIncompleteCompl( a, x ) = p.
280 real gammaIncompleteComplInverse(real a
, real p
)
282 assert(p
>= 0 && p
<= 1);
286 return std
.internal
.math
.gammafunction
.gammaIncompleteComplInv(a
, p
);
290 /* ***********************************************
291 * ERROR FUNCTIONS & NORMAL DISTRIBUTION *
292 * ***********************************************/
298 * erf(x) = 2/ $(SQRT)($(PI))
299 * $(INTEGRATE 0, x) exp( - $(POWER t, 2)) dt
301 * The magnitude of x is limited to about 106.56 for IEEE 80-bit
302 * arithmetic; 1 or -1 is returned outside this range.
306 return std
.internal
.math
.errorfunction
.erf(x
);
309 /** Complementary error function
311 * erfc(x) = 1 - erf(x)
312 * = 2/ $(SQRT)($(PI))
313 * $(INTEGRATE x, $(INFIN)) exp( - $(POWER t, 2)) dt
315 * This function has high relative accuracy for
316 * values of x far from zero. (For values near zero, use erf(x)).
320 return std
.internal
.math
.errorfunction
.erfc(x
);
324 /** Normal distribution function.
326 * The normal (or Gaussian, or bell-shaped) distribution is
329 * normalDist(x) = 1/$(SQRT)(2$(PI)) $(INTEGRATE -$(INFIN), x) exp( - $(POWER t, 2)/2) dt
330 * = 0.5 + 0.5 * erf(x/sqrt(2))
331 * = 0.5 * erfc(- x/sqrt(2))
333 * To maintain accuracy at values of x near 1.0, use
334 * normalDistribution(x) = 1.0 - normalDistribution(-x).
337 * $(LINK http://www.netlib.org/cephes/ldoubdoc.html),
338 * G. Marsaglia, "Evaluating the Normal Distribution",
339 * Journal of Statistical Software <b>11</b>, (July 2004).
341 real normalDistribution(real x
)
343 return std
.internal
.math
.errorfunction
.normalDistributionImpl(x
);
346 /** Inverse of Normal distribution function
348 * Returns the argument, x, for which the area under the
349 * Normal probability density function (integrated from
350 * minus infinity to x) is equal to p.
352 * Note: This function is only implemented to 80 bit precision.
354 real normalDistributionInverse(real p
)
356 assert(p
>= 0.0L && p
<= 1.0L, "Domain error");
360 return std
.internal
.math
.errorfunction
.normalDistributionInvImpl(p
);