std.mathspecial source code

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  *      PSI = &Psi;
33  *      THETA = &theta;
34  *      INTEGRAL = &#8747;
35  *      INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
36  *      POWER = $1<sup>$2</sup>
37  *      SUB = $1<sub>$2</sub>
38  *      BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
39  *      CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
40  *      CEIL = &#8968;$1&#8969;
41  *      PLUSMN = &plusmn;
42  *      MNPLUS = &mnplus;
43  *      INFIN = &infin;
44  *      PLUSMNINF = &plusmn;&infin;
45  *      MNPLUSINF = &mnplus;&infin;
46  *      PI = &pi;
47  *      LT = &lt;
48  *      LE = &le;
49  *      GT = &gt;
50  *      SQRT = &radic;
51  *      HALF = &frac12;
52  *      COMPLEX = &#8450;
53  *
54  * Copyright: Based on the CEPHES math library, which is
55  *            Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
56  * License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
57  * Authors:   Stephen L. Moshier (original C code). Conversion to D by Don Clugston
58  * Source:    $(PHOBOSSRC std/mathspecial.d)
59  */
60 module std.mathspecial;
61 import std.internal.math.errorfunction;
62 import std.internal.math.gammafunction;
63 public import std.math;
64 
65 /* ***********************************************
66  *            GAMMA AND RELATED FUNCTIONS        *
67  * ***********************************************/
68 
69 pure:
70 nothrow:
71 @safe:
72 @nogc:
73 
74 /** The Gamma function, $(GAMMA)(x)
75  *
76  *  $(GAMMA)(x) is a generalisation of the factorial function
77  *  to real and complex numbers.
78  *  Like x!, $(GAMMA)(x+1) = x * $(GAMMA)(x).
79  *
80  *  Mathematically, if z.re > 0 then
81  *   $(GAMMA)(z) = $(INTEGRATE 0, $(INFIN)) $(POWER t, z-1)$(POWER e, -t) dt
82  *
83  *  $(TABLE_SV
84  *    $(SVH  x,           $(GAMMA)(x) )
85  *    $(SV  $(NAN),       $(NAN)      )
86  *    $(SV  $(PLUSMN)0.0, $(PLUSMNINF))
87  *    $(SV integer > 0,   (x-1)!      )
88  *    $(SV integer < 0,   $(NAN)      )
89  *    $(SV +$(INFIN),      +$(INFIN)   )
90  *    $(SV -$(INFIN),      $(NAN)      )
91  *  )
92  */
93 real gamma(real x)
94 {
95     return std.internal.math.gammafunction.gamma(x);
96 }
97 
98 /** Natural logarithm of the gamma function, $(GAMMA)(x)
99  *
100  * Returns the base e (2.718...) logarithm of the absolute
101  * value of the gamma function of the argument.
102  *
103  * For reals, logGamma is equivalent to log(fabs(gamma(x))).
104  *
105  *  $(TABLE_SV
106  *    $(SVH  x,             logGamma(x)   )
107  *    $(SV  $(NAN),         $(NAN)      )
108  *    $(SV integer <= 0,    +$(INFIN)    )
109  *    $(SV $(PLUSMNINF),    +$(INFIN)    )
110  *  )
111  */
112 real logGamma(real x)
113 {
114     return std.internal.math.gammafunction.logGamma(x);
115 }
116 
117 /** The sign of $(GAMMA)(x).
118  *
119  * Params:
120  *   x = the argument of $(GAMMA)
121  *
122  * Returns:
123  *   -1 if $(GAMMA)(x) < 0, +1 if $(GAMMA)(x) > 0, and $(NAN) if $(GAMMA)(x)
124  *   does not exist.
125  *
126  * Note:
127  *   This function can be used in conjunction with `logGamma` to evaluate
128  *   $(GAMMA)(x) when `gamma(x)` is too large to be represented as a `real`.
129  */
130 pragma(inline, true) real sgnGamma(real x)
131 {
132     return std.internal.math.gammafunction.sgnGamma(x);
133 }
134 
135 ///
136 @safe unittest
137 {
138     assert(sgnGamma(10_000) == 1);
139 }
140 
141 /** Beta function, B(x,y)
142  *
143  * Mathematically, if x $(GT) 0 and y $(GT) 0 then
144  * B(x,y) = $(INTEGRATE 0, 1)$(POWER t, x-1)$(POWER (l-t), y-1)dt. Through analytic continuation, it
145  * is extended to $(COMPLEX)$(SUP 2) where it can be expressed in terms of $(GAMMA)(z).
146  *
147  * B(x,y) = $(GAMMA)(x)$(GAMMA)(y) / $(GAMMA)(x+y).
148  *
149  * This implementation restricts x and y to the set of real numbers.
150  *
151  * Params:
152  *   x = the first argument of B
153  *   y = the second argument of B
154  *
155  * Returns:
156  *   It returns B(x,y) if it can be computed, otherwise $(NAN).
157  *
158  * $(TABLE_SV
159  *   $(TR $(TH x)                                   $(TH y)                $(TH beta(x, y))   )
160  *   $(TR $(TD $(NAN))                              $(TD y)                $(TD $(NAN))       )
161  *   $(TR $(TD -$(INFIN))                           $(TD y)                $(TD $(NAN))       )
162  *   $(TR $(TD integer $(LT) 0)                     $(TD y)                $(TD $(NAN))       )
163  *   $(TR $(TD noninteger and x+y even $(LE) 0)     $(TD noninteger)       $(TD -0)           )
164  *   $(TR $(TD noninteger and x+y odd $(LE) 0)      $(TD noninteger)       $(TD +0)           )
165  *   $(TR $(TD +0)                                  $(TD positive finite)  $(TD +$(INFIN))    )
166  *   $(TR $(TD +0)                                  $(TD +$(INFIN))        $(TD $(NAN))       )
167  *   $(TR $(TD $(GT) 0)                             $(TD +$(INFIN))        $(TD +0)           )
168  *   $(TR $(TD -0)                                  $(TD +0)               $(TD $(NAN))       )
169  *   $(TR $(TD -0)                                  $(TD $(GT) 0)          $(TD -$(INFIN))    )
170  *   $(TR $(TD noninteger $(LT) 0, $(CEIL x) odd)   $(TD +$(INFIN))        $(TD -$(INFIN))    )
171  *   $(TR $(TD noninteger $(LT) 0, $(CEIL x) even)  $(TD +$(INFIN))        $(TD +$(INFIN))    )
172  *   $(TR $(TD noninteger $(LT) 0)                  $(TD $(PLUSMN)0)       $(TD $(PLUSMNINF)) )
173  * )
174  *
175  * Since B(x,y) = B(y,x), if the table states that beta(x, y) is a special value, then beta(y, x) is
176  * one as well.
177  */
178 pragma(inline, true) real beta(real x, real y)
179 {
180     return std.internal.math.gammafunction.beta(x, y);
181 }
182 
183 ///
184 @safe unittest
185 {
186     assert(beta(1, 2) == 0.5);
187     assert(isIdentical(beta(NaN(0xABC), 4), NaN(0xABC)));
188     assert(beta(3, 4) == beta(4, 3));
189     assert(isNaN(beta(-real.infinity, +0.)));
190     assert(isNaN(beta(-1, 2)));
191     assert(beta(-0.5, 0.5) is -0.0L);
192     assert(beta(-1.5, 0.5) is +0.0L);
193     assert(beta(+0., +0.) == +real.infinity);
194     assert(isNaN(beta(+0., +real.infinity)));
195     assert(beta(1, +real.infinity) is +0.0L);
196     assert(isNaN(beta(-0., +0.)));
197     assert(beta(-0., nextUp(+0.0L)) == -real.infinity);
198     assert(beta(-0.5, +real.infinity) == -real.infinity);
199     assert(beta(nextDown(-1.0L), real.infinity) == real.infinity);
200     assert(beta(nextDown(-0.0L), +0.) == +real.infinity);
201     assert(beta(-0.5, -0.) == -real.infinity);
202 }
203 
204 /** Digamma function, $(PSI)(x)
205  *
206  * $(PSI)(x), is the logarithmic derivative of the gamma function, $(GAMMA)(x).
207  *
208  * $(PSI)(x) = $(SUP d)$(SUB /, dx) ln|$(GAMMA)(x)|  (the derivative of `logGamma(x)`)
209  *
210  * Params:
211  *   x = the domain value
212  *
213  * Returns:
214  *   It returns $(PSI)(x).
215  *
216  * $(TABLE_SV
217  *   $(SVH x,             digamma(x)   )
218  *   $(SV  integer < 0,   $(NAN)       )
219  *   $(SV  $(PLUSMN)0.0,  $(MNPLUSINF) )
220  *   $(SV  +$(INFIN),     +$(INFIN)    )
221  *   $(SV  -$(INFIN),     $(NAN)       )
222  *   $(SV  $(NAN),        $(NAN)       )
223  * )
224  *
225  * See_Also: $(LREF logmdigamma), $(LREF logmdigammaInverse).
226  */
227 real digamma(real x)
228 {
229     return std.internal.math.gammafunction.digamma(x);
230 }
231 
232 ///
233 @safe unittest
234 {
235     const euler = 0.57721_56649_01532_86060_65121L;
236 
237     assert(isClose(digamma(1), -euler));
238     assert(digamma(+0.) == -real.infinity);
239     assert(digamma(-0.) == +real.infinity);
240     assert(digamma(+real.infinity) == +real.infinity);
241     assert(isNaN(digamma(-1)));
242     assert(isNaN(digamma(-real.infinity)));
243 }
244 
245 /** Log Minus Digamma function
246  *
247  *  logmdigamma(x) = log(x) - digamma(x)
248  *
249  *  See_Also: $(LREF digamma), $(LREF logmdigammaInverse).
250  */
251 real logmdigamma(real x)
252 {
253     return std.internal.math.gammafunction.logmdigamma(x);
254 }
255 
256 /** Inverse of the Log Minus Digamma function
257  *
258  *  Given y, the function finds x such log(x) - digamma(x) = y.
259  *
260  *  See_Also: $(LREF logmdigamma), $(LREF digamma).
261  */
262 real logmdigammaInverse(real x)
263 {
264     return std.internal.math.gammafunction.logmdigammaInverse(x);
265 }
266 
267 /** Regularized incomplete beta function $(SUB I, x)(a,b)
268  *
269  * Mathematically, if a and b are positive real numbers, and 0 $(LE) x $(LE) 1, then
270  * $(SUB I, x)(a,b) = $(INTEGRATE 0, x)$(POWER t, a-1)$(POWER (1-t), b-1)dt/B(a,b) where B is the
271  * beta function. It is also the cumulative distribution function of the beta distribution.
272  *
273  * `betaIncomplete(a, b, x)` evaluates $(SUB I, `x`)(`a`,`b`).
274  *
275  * Params:
276  *   a = the first argument of B, must be positive
277  *   b = the second argument of B, must be positive
278  *   x = the fraction of integration completion from below, 0 $(LE) x $(LE) 1
279  *
280  * Returns:
281  *   It returns $(SUB I, x)(a,b), an element of [0,1].
282  *
283  * Note:
284  *   The integral is evaluated by a continued fraction expansion or, when `b * x` is small, by a
285  *   power series.
286  *
287  * See_Also: $(LREF beta) $(LREF betaIncompleteCompl)
288  */
289 real betaIncomplete(real a, real b, real x )
290 {
291     return std.internal.math.gammafunction.betaIncomplete(a, b, x);
292 }
293 
294 /** Regularized incomplete beta function complement $(SUB I$(SUP C), x)(a,b)
295  *
296  * Mathematically, if a $(GT) 0, b $(GT) 0, and 0 $(LE) x $(LE) 1, then
297  * $(SUB I$(SUP C), x)(a,b) = $(INTEGRATE x, 1)$(POWER t, a-1)$(POWER (1-t), b-1)dt/B(a,b) where B
298  * is the beta function. It is also the complement of the cumulative distribution function of the
299  * beta distribution. It can be shown that $(SUB I$(SUP C), x)(a,b) = $(SUB I, 1-x)(b,a).
300  *
301  * `betaIncompleteCompl(a, b, x)` evaluates $(SUB I$(SUP C), `x`)(`a`,`b`).
302  *
303  * Params:
304  *   a = the first argument of B, must be positive
305  *   b = the second argument of B, must be positive
306  *   x = the fraction of integration completion from above, 0 $(LE) x $(LE) 1
307  *
308  * Returns:
309  *   It returns $(SUB I$(SUP C), x)(a,b), an element of [0,1].
310  *
311  * See_Also: $(LREF beta) $(LREF betaIncomplete)
312  */
313 real betaIncompleteCompl(real a, real b, real x)
314 in
315 {
316     // allow NaN input to pass through so that it can be addressed by the
317     // internal NaN payload propagation logic
318     if (!isNaN(a) && !isNaN(b) && !isNaN(x))
319     {
320         assert(signbit(a) == 0, "a must be positive");
321         assert(signbit(b) == 0, "b must be positive");
322         assert(x >= 0 && x <= 1, "x must be in [0, 1]");
323     }
324 }
325 do
326 {
327     return std.internal.math.gammafunction.betaIncomplete(b, a, 1-x);
328 }
329 
330 ///
331 @safe unittest
332 {
333     assert(betaIncompleteCompl(.1, .2, 0) == betaIncomplete(.2, .1, 1));
334 }
335 
336 /** Inverse of incomplete beta integral
337  *
338  * Given y, the function finds x such that
339  *
340  *  betaIncomplete(a, b, x) == y
341  *
342  *  Newton iterations or interval halving is used.
343  */
344 real betaIncompleteInverse(real a, real b, real y )
345 {
346     return std.internal.math.gammafunction.betaIncompleteInv(a, b, y);
347 }
348 
349 /** Incomplete gamma integral and its complement
350  *
351  * These functions are defined by
352  *
353  *   gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
354  *
355  *  gammaIncompleteCompl(a,x)   =   1 - gammaIncomplete(a,x)
356  * = ($(INTEGRATE x, $(INFIN)) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
357  *
358  * In this implementation both arguments must be positive.
359  * The integral is evaluated by either a power series or
360  * continued fraction expansion, depending on the relative
361  * values of a and x.
362  */
363 real gammaIncomplete(real a, real x )
364 in
365 {
366    assert(x >= 0);
367    assert(a > 0);
368 }
369 do
370 {
371     return std.internal.math.gammafunction.gammaIncomplete(a, x);
372 }
373 
374 /** ditto */
375 real gammaIncompleteCompl(real a, real x )
376 in
377 {
378    assert(x >= 0);
379    assert(a > 0);
380 }
381 do
382 {
383     return std.internal.math.gammafunction.gammaIncompleteCompl(a, x);
384 }
385 
386 /** Inverse of complemented incomplete gamma integral
387  *
388  * Given a and p, the function finds x such that
389  *
390  *  gammaIncompleteCompl( a, x ) = p.
391  */
392 real gammaIncompleteComplInverse(real a, real p)
393 in
394 {
395   assert(p >= 0 && p <= 1);
396   assert(a > 0);
397 }
398 do
399 {
400     return std.internal.math.gammafunction.gammaIncompleteComplInv(a, p);
401 }
402 
403 
404 /* ***********************************************
405  *     ERROR FUNCTIONS & NORMAL DISTRIBUTION     *
406  * ***********************************************/
407 
408  /** Error function
409  *
410  * The integral is
411  *
412  *  erf(x) =  2/ $(SQRT)($(PI))
413  *     $(INTEGRATE 0, x) exp( - $(POWER t, 2)) dt
414  *
415  * The magnitude of x is limited to about 106.56 for IEEE 80-bit
416  * arithmetic; 1 or -1 is returned outside this range.
417  */
418 real erf(real x)
419 {
420     return std.internal.math.errorfunction.erf(x);
421 }
422 
423 /** Complementary error function
424  *
425  * erfc(x) = 1 - erf(x)
426  *         = 2/ $(SQRT)($(PI))
427  *     $(INTEGRATE x, $(INFIN)) exp( - $(POWER t, 2)) dt
428  *
429  * This function has high relative accuracy for
430  * values of x far from zero. (For values near zero, use erf(x)).
431  */
432 real erfc(real x)
433 {
434     return std.internal.math.errorfunction.erfc(x);
435 }
436 
437 
438 /** Standard normal distribution function.
439  *
440  * The normal (or Gaussian, or bell-shaped) distribution is
441  * defined as:
442  *
443  * normalDist(x) = 1/$(SQRT)(2$(PI)) $(INTEGRATE -$(INFIN), x) exp( - $(POWER t, 2)/2) dt
444  *   = 0.5 + 0.5 * erf(x/sqrt(2))
445  *   = 0.5 * erfc(- x/sqrt(2))
446  *
447  * To maintain accuracy at values of x near 1.0, use
448  *      normalDistribution(x) = 1.0 - normalDistribution(-x).
449  *
450  * References:
451  * $(LINK http://www.netlib.org/cephes/ldoubdoc.html),
452  * G. Marsaglia, "Evaluating the Normal Distribution",
453  * Journal of Statistical Software <b>11</b>, (July 2004).
454  */
455 real normalDistribution(real x)
456 {
457     return std.internal.math.errorfunction.normalDistributionImpl(x);
458 }
459 
460 /** Inverse of Standard normal distribution function
461  *
462  * Returns the argument, x, for which the area under the
463  * Normal probability density function (integrated from
464  * minus infinity to x) is equal to p.
465  *
466  * Note: This function is only implemented to 80 bit precision.
467  */
468 real normalDistributionInverse(real p)
469 in
470 {
471   assert(p >= 0.0L && p <= 1.0L, "Domain error");
472 }
473 do
474 {
475     return std.internal.math.errorfunction.normalDistributionInvImpl(p);
476 }