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 */
55module std.mathspecial;
56import std.internal.math.errorfunction;
57import std.internal.math.gammafunction;
58public import std.math;
59
60/* ***********************************************
61 *            GAMMA AND RELATED FUNCTIONS        *
62 * ***********************************************/
63
64pure:
65nothrow:
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 */
88real 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 */
107real 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 */
120real 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 */
155real 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 */
177real 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 */
188real 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 */
199real 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 */
225real 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 */
238real 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 */
257real gammaIncomplete(real a, real x )
258in
259{
260   assert(x >= 0);
261   assert(a > 0);
262}
263do
264{
265    return std.internal.math.gammafunction.gammaIncomplete(a, x);
266}
267
268/** ditto */
269real gammaIncompleteCompl(real a, real x )
270in
271{
272   assert(x >= 0);
273   assert(a > 0);
274}
275do
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 */
286real gammaIncompleteComplInverse(real a, real p)
287in
288{
289  assert(p >= 0 && p <= 1);
290  assert(a > 0);
291}
292do
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 */
312real 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 */
326real 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 */
349real 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 */
362real normalDistributionInverse(real p)
363in
364{
365  assert(p >= 0.0L && p <= 1.0L, "Domain error");
366}
367do
368{
369    return std.internal.math.errorfunction.normalDistributionInvImpl(p);
370}
371