The OpenD Programming Language

1 // Written in the D programming language.
2 
3 /**
4  * Mathematical Special Functions
5  *
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.
9  *
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.
16  *
17  * Status:
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.
21  *
22  * Macros:
23  *      TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
24  *              <caption>Special Values</caption>
25  *              $0</table>
26  *      SVH = $(TR $(TH $1) $(TH $2))
27  *      SV  = $(TR $(TD $1) $(TD $2))
28  *
29  *      NAN = $(RED NAN)
30  *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
31  *      GAMMA = &#915;
32  *      THETA = &theta;
33  *      INTEGRAL = &#8747;
34  *      INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
35  *      POWER = $1<sup>$2</sup>
36  *      SUB = $1<sub>$2</sub>
37  *      BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
38  *      CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
39  *      PLUSMN = &plusmn;
40  *      INFIN = &infin;
41  *      PLUSMNINF = &plusmn;&infin;
42  *      PI = &pi;
43  *      LT = &lt;
44  *      GT = &gt;
45  *      SQRT = &radic;
46  *      HALF = &frac12;
47  *
48  *
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)
54  */
55 module std.mathspecial;
56 import std.internal.math.errorfunction;
57 import std.internal.math.gammafunction;
58 public import std.math;
59 
60 /* ***********************************************
61  *            GAMMA AND RELATED FUNCTIONS        *
62  * ***********************************************/
63 
64 pure:
65 nothrow:
66 @safe:
67 @nogc:
68 
69 /** The Gamma function, $(GAMMA)(x)
70  *
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).
74  *
75  *  Mathematically, if z.re > 0 then
76  *   $(GAMMA)(z) = $(INTEGRATE 0, $(INFIN)) $(POWER t, z-1)$(POWER e, -t) dt
77  *
78  *  $(TABLE_SV
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)      )
86  *  )
87  */
88 real gamma(real x)
89 {
90     return std.internal.math.gammafunction.gamma(x);
91 }
92 
93 /** Natural logarithm of the gamma function, $(GAMMA)(x)
94  *
95  * Returns the base e (2.718...) logarithm of the absolute
96  * value of the gamma function of the argument.
97  *
98  * For reals, logGamma is equivalent to log(fabs(gamma(x))).
99  *
100  *  $(TABLE_SV
101  *    $(SVH  x,             logGamma(x)   )
102  *    $(SV  $(NAN),         $(NAN)      )
103  *    $(SV integer <= 0,    +$(INFIN)    )
104  *    $(SV $(PLUSMNINF),    +$(INFIN)    )
105  *  )
106  */
107 real logGamma(real x)
108 {
109     return std.internal.math.gammafunction.logGamma(x);
110 }
111 
112 /** The sign of $(GAMMA)(x).
113  *
114  * Returns -1 if $(GAMMA)(x) < 0,  +1 if $(GAMMA)(x) > 0,
115  * $(NAN) if sign is indeterminate.
116  *
117  * Note that this function can be used in conjunction with logGamma(x) to
118  * evaluate gamma for very large values of x.
119  */
120 real sgnGamma(real x)
121 {
122     import core.math : rndtol;
123     /* Author: Don Clugston. */
124     if (isNaN(x)) return x;
125     if (x > 0) return 1.0;
126     if (x < -1/real.epsilon)
127     {
128         // Large negatives lose all precision
129         return real.nan;
130     }
131     long n = rndtol(x);
132     if (x == n)
133     {
134         return x == 0 ?  copysign(1, x) : real.nan;
135     }
136     return n & 1 ? 1.0 : -1.0;
137 }
138 
139 @safe unittest
140 {
141     assert(sgnGamma(5.0) == 1.0);
142     assert(isNaN(sgnGamma(-3.0)));
143     assert(sgnGamma(-0.1) == -1.0);
144     assert(sgnGamma(-55.1) == 1.0);
145     assert(isNaN(sgnGamma(-real.infinity)));
146     assert(isIdentical(sgnGamma(NaN(0xABC)), NaN(0xABC)));
147 }
148 
149 /** Beta function
150  *
151  * The beta function is defined as
152  *
153  * beta(x, y) = ($(GAMMA)(x) * $(GAMMA)(y)) / $(GAMMA)(x + y)
154  */
155 real beta(real x, real y)
156 {
157     if ((x+y)> MAXGAMMA)
158     {
159         return exp(logGamma(x) + logGamma(y) - logGamma(x+y));
160     } else return gamma(x) * gamma(y) / gamma(x+y);
161 }
162 
163 @safe unittest
164 {
165     assert(isIdentical(beta(NaN(0xABC), 4), NaN(0xABC)));
166     assert(isIdentical(beta(2, NaN(0xABC)), NaN(0xABC)));
167 }
168 
169 /** Digamma function
170  *
171  *  The digamma function is the logarithmic derivative of the gamma function.
172  *
173  *  digamma(x) = d/dx logGamma(x)
174  *
175  *  See_Also: $(LREF logmdigamma), $(LREF logmdigammaInverse).
176  */
177 real digamma(real x)
178 {
179     return std.internal.math.gammafunction.digamma(x);
180 }
181 
182 /** Log Minus Digamma function
183  *
184  *  logmdigamma(x) = log(x) - digamma(x)
185  *
186  *  See_Also: $(LREF digamma), $(LREF logmdigammaInverse).
187  */
188 real logmdigamma(real x)
189 {
190     return std.internal.math.gammafunction.logmdigamma(x);
191 }
192 
193 /** Inverse of the Log Minus Digamma function
194  *
195  *  Given y, the function finds x such log(x) - digamma(x) = y.
196  *
197  *  See_Also: $(LREF logmdigamma), $(LREF digamma).
198  */
199 real logmdigammaInverse(real x)
200 {
201     return std.internal.math.gammafunction.logmdigammaInverse(x);
202 }
203 
204 /** Incomplete beta integral
205  *
206  * Returns regularized incomplete beta integral of the arguments, evaluated
207  * from zero to x. The regularized incomplete beta function is defined as
208  *
209  * betaIncomplete(a, b, x) = $(GAMMA)(a + b) / ( $(GAMMA)(a) $(GAMMA)(b) ) *
210  * $(INTEGRATE 0, x) $(POWER t, a-1)$(POWER (1-t), b-1) dt
211  *
212  * and is the same as the cumulative distribution function of the Beta
213  * distribution.
214  *
215  * The domain of definition is 0 <= x <= 1.  In this
216  * implementation a and b are restricted to positive values.
217  * The integral from x to 1 may be obtained by the symmetry
218  * relation
219  *
220  *    betaIncompleteCompl(a, b, x )  =  betaIncomplete( b, a, 1-x )
221  *
222  * The integral is evaluated by a continued fraction expansion
223  * or, when b * x is small, by a power series.
224  */
225 real betaIncomplete(real a, real b, real x )
226 {
227     return std.internal.math.gammafunction.betaIncomplete(a, b, x);
228 }
229 
230 /** Inverse of incomplete beta integral
231  *
232  * Given y, the function finds x such that
233  *
234  *  betaIncomplete(a, b, x) == y
235  *
236  *  Newton iterations or interval halving is used.
237  */
238 real betaIncompleteInverse(real a, real b, real y )
239 {
240     return std.internal.math.gammafunction.betaIncompleteInv(a, b, y);
241 }
242 
243 /** Incomplete gamma integral and its complement
244  *
245  * These functions are defined by
246  *
247  *   gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
248  *
249  *  gammaIncompleteCompl(a,x)   =   1 - gammaIncomplete(a,x)
250  * = ($(INTEGRATE x, $(INFIN)) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
251  *
252  * In this implementation both arguments must be positive.
253  * The integral is evaluated by either a power series or
254  * continued fraction expansion, depending on the relative
255  * values of a and x.
256  */
257 real gammaIncomplete(real a, real x )
258 in
259 {
260    assert(x >= 0);
261    assert(a > 0);
262 }
263 do
264 {
265     return std.internal.math.gammafunction.gammaIncomplete(a, x);
266 }
267 
268 /** ditto */
269 real gammaIncompleteCompl(real a, real x )
270 in
271 {
272    assert(x >= 0);
273    assert(a > 0);
274 }
275 do
276 {
277     return std.internal.math.gammafunction.gammaIncompleteCompl(a, x);
278 }
279 
280 /** Inverse of complemented incomplete gamma integral
281  *
282  * Given a and p, the function finds x such that
283  *
284  *  gammaIncompleteCompl( a, x ) = p.
285  */
286 real gammaIncompleteComplInverse(real a, real p)
287 in
288 {
289   assert(p >= 0 && p <= 1);
290   assert(a > 0);
291 }
292 do
293 {
294     return std.internal.math.gammafunction.gammaIncompleteComplInv(a, p);
295 }
296 
297 
298 /* ***********************************************
299  *     ERROR FUNCTIONS & NORMAL DISTRIBUTION     *
300  * ***********************************************/
301 
302  /** Error function
303  *
304  * The integral is
305  *
306  *  erf(x) =  2/ $(SQRT)($(PI))
307  *     $(INTEGRATE 0, x) exp( - $(POWER t, 2)) dt
308  *
309  * The magnitude of x is limited to about 106.56 for IEEE 80-bit
310  * arithmetic; 1 or -1 is returned outside this range.
311  */
312 real erf(real x)
313 {
314     return std.internal.math.errorfunction.erf(x);
315 }
316 
317 /** Complementary error function
318  *
319  * erfc(x) = 1 - erf(x)
320  *         = 2/ $(SQRT)($(PI))
321  *     $(INTEGRATE x, $(INFIN)) exp( - $(POWER t, 2)) dt
322  *
323  * This function has high relative accuracy for
324  * values of x far from zero. (For values near zero, use erf(x)).
325  */
326 real erfc(real x)
327 {
328     return std.internal.math.errorfunction.erfc(x);
329 }
330 
331 
332 /** Standard normal distribution function.
333  *
334  * The normal (or Gaussian, or bell-shaped) distribution is
335  * defined as:
336  *
337  * normalDist(x) = 1/$(SQRT)(2$(PI)) $(INTEGRATE -$(INFIN), x) exp( - $(POWER t, 2)/2) dt
338  *   = 0.5 + 0.5 * erf(x/sqrt(2))
339  *   = 0.5 * erfc(- x/sqrt(2))
340  *
341  * To maintain accuracy at values of x near 1.0, use
342  *      normalDistribution(x) = 1.0 - normalDistribution(-x).
343  *
344  * References:
345  * $(LINK http://www.netlib.org/cephes/ldoubdoc.html),
346  * G. Marsaglia, "Evaluating the Normal Distribution",
347  * Journal of Statistical Software <b>11</b>, (July 2004).
348  */
349 real normalDistribution(real x)
350 {
351     return std.internal.math.errorfunction.normalDistributionImpl(x);
352 }
353 
354 /** Inverse of Standard normal distribution function
355  *
356  * Returns the argument, x, for which the area under the
357  * Normal probability density function (integrated from
358  * minus infinity to x) is equal to p.
359  *
360  * Note: This function is only implemented to 80 bit precision.
361  */
362 real normalDistributionInverse(real p)
363 in
364 {
365   assert(p >= 0.0L && p <= 1.0L, "Domain error");
366 }
367 do
368 {
369     return std.internal.math.errorfunction.normalDistributionInvImpl(p);
370 }