|
145
|
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 = Γ
|
|
|
32 * THETA = θ
|
|
|
33 * INTEGRAL = ∫
|
|
|
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 ))
|
|
|
39 * PLUSMN = ±
|
|
|
40 * INFIN = ∞
|
|
|
41 * PLUSMNINF = ±∞
|
|
|
42 * PI = π
|
|
|
43 * LT = <
|
|
|
44 * GT = >
|
|
|
45 * SQRT = √
|
|
|
46 * HALF = ½
|
|
|
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 /* Author: Don Clugston. */
|
|
|
123 if (isNaN(x)) return x;
|
|
|
124 if (x > 0) return 1.0;
|
|
|
125 if (x < -1/real.epsilon)
|
|
|
126 {
|
|
|
127 // Large negatives lose all precision
|
|
|
128 return real.nan;
|
|
|
129 }
|
|
|
130 long n = rndtol(x);
|
|
|
131 if (x == n)
|
|
|
132 {
|
|
|
133 return x == 0 ? copysign(1, x) : real.nan;
|
|
|
134 }
|
|
|
135 return n & 1 ? 1.0 : -1.0;
|
|
|
136 }
|
|
|
137
|
|
|
138 @safe unittest
|
|
|
139 {
|
|
|
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)));
|
|
|
146 }
|
|
|
147
|
|
|
148 /** Beta function
|
|
|
149 *
|
|
|
150 * The beta function is defined as
|
|
|
151 *
|
|
|
152 * beta(x, y) = ($(GAMMA)(x) * $(GAMMA)(y)) / $(GAMMA)(x + y)
|
|
|
153 */
|
|
|
154 real beta(real x, real y)
|
|
|
155 {
|
|
|
156 if ((x+y)> MAXGAMMA)
|
|
|
157 {
|
|
|
158 return exp(logGamma(x) + logGamma(y) - logGamma(x+y));
|
|
|
159 } else return gamma(x) * gamma(y) / gamma(x+y);
|
|
|
160 }
|
|
|
161
|
|
|
162 @safe unittest
|
|
|
163 {
|
|
|
164 assert(isIdentical(beta(NaN(0xABC), 4), NaN(0xABC)));
|
|
|
165 assert(isIdentical(beta(2, NaN(0xABC)), NaN(0xABC)));
|
|
|
166 }
|
|
|
167
|
|
|
168 /** Digamma function
|
|
|
169 *
|
|
|
170 * The digamma function is the logarithmic derivative of the gamma function.
|
|
|
171 *
|
|
|
172 * digamma(x) = d/dx logGamma(x)
|
|
|
173 *
|
|
|
174 * See_Also: $(LREF logmdigamma), $(LREF logmdigammaInverse).
|
|
|
175 */
|
|
|
176 real digamma(real x)
|
|
|
177 {
|
|
|
178 return std.internal.math.gammafunction.digamma(x);
|
|
|
179 }
|
|
|
180
|
|
|
181 /** Log Minus Digamma function
|
|
|
182 *
|
|
|
183 * logmdigamma(x) = log(x) - digamma(x)
|
|
|
184 *
|
|
|
185 * See_Also: $(LREF digamma), $(LREF logmdigammaInverse).
|
|
|
186 */
|
|
|
187 real logmdigamma(real x)
|
|
|
188 {
|
|
|
189 return std.internal.math.gammafunction.logmdigamma(x);
|
|
|
190 }
|
|
|
191
|
|
|
192 /** Inverse of the Log Minus Digamma function
|
|
|
193 *
|
|
|
194 * Given y, the function finds x such log(x) - digamma(x) = y.
|
|
|
195 *
|
|
|
196 * See_Also: $(LREF logmdigamma), $(LREF digamma).
|
|
|
197 */
|
|
|
198 real logmdigammaInverse(real x)
|
|
|
199 {
|
|
|
200 return std.internal.math.gammafunction.logmdigammaInverse(x);
|
|
|
201 }
|
|
|
202
|
|
|
203 /** Incomplete beta integral
|
|
|
204 *
|
|
|
205 * Returns incomplete beta integral of the arguments, evaluated
|
|
|
206 * from zero to x. The regularized incomplete beta function is defined as
|
|
|
207 *
|
|
|
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
|
|
|
210 *
|
|
|
211 * and is the same as the the cumulative distribution function.
|
|
|
212 *
|
|
|
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
|
|
|
216 * relation
|
|
|
217 *
|
|
|
218 * betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x )
|
|
|
219 *
|
|
|
220 * The integral is evaluated by a continued fraction expansion
|
|
|
221 * or, when b * x is small, by a power series.
|
|
|
222 */
|
|
|
223 real betaIncomplete(real a, real b, real x )
|
|
|
224 {
|
|
|
225 return std.internal.math.gammafunction.betaIncomplete(a, b, x);
|
|
|
226 }
|
|
|
227
|
|
|
228 /** Inverse of incomplete beta integral
|
|
|
229 *
|
|
|
230 * Given y, the function finds x such that
|
|
|
231 *
|
|
|
232 * betaIncomplete(a, b, x) == y
|
|
|
233 *
|
|
|
234 * Newton iterations or interval halving is used.
|
|
|
235 */
|
|
|
236 real betaIncompleteInverse(real a, real b, real y )
|
|
|
237 {
|
|
|
238 return std.internal.math.gammafunction.betaIncompleteInv(a, b, y);
|
|
|
239 }
|
|
|
240
|
|
|
241 /** Incomplete gamma integral and its complement
|
|
|
242 *
|
|
|
243 * These functions are defined by
|
|
|
244 *
|
|
|
245 * gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
|
|
|
246 *
|
|
|
247 * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x)
|
|
|
248 * = ($(INTEGRATE x, $(INFIN)) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
|
|
|
249 *
|
|
|
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
|
|
|
253 * values of a and x.
|
|
|
254 */
|
|
|
255 real gammaIncomplete(real a, real x )
|
|
|
256 in {
|
|
|
257 assert(x >= 0);
|
|
|
258 assert(a > 0);
|
|
|
259 }
|
|
|
260 body {
|
|
|
261 return std.internal.math.gammafunction.gammaIncomplete(a, x);
|
|
|
262 }
|
|
|
263
|
|
|
264 /** ditto */
|
|
|
265 real gammaIncompleteCompl(real a, real x )
|
|
|
266 in {
|
|
|
267 assert(x >= 0);
|
|
|
268 assert(a > 0);
|
|
|
269 }
|
|
|
270 body {
|
|
|
271 return std.internal.math.gammafunction.gammaIncompleteCompl(a, x);
|
|
|
272 }
|
|
|
273
|
|
|
274 /** Inverse of complemented incomplete gamma integral
|
|
|
275 *
|
|
|
276 * Given a and p, the function finds x such that
|
|
|
277 *
|
|
|
278 * gammaIncompleteCompl( a, x ) = p.
|
|
|
279 */
|
|
|
280 real gammaIncompleteComplInverse(real a, real p)
|
|
|
281 in {
|
|
|
282 assert(p >= 0 && p <= 1);
|
|
|
283 assert(a > 0);
|
|
|
284 }
|
|
|
285 body {
|
|
|
286 return std.internal.math.gammafunction.gammaIncompleteComplInv(a, p);
|
|
|
287 }
|
|
|
288
|
|
|
289
|
|
|
290 /* ***********************************************
|
|
|
291 * ERROR FUNCTIONS & NORMAL DISTRIBUTION *
|
|
|
292 * ***********************************************/
|
|
|
293
|
|
|
294 /** Error function
|
|
|
295 *
|
|
|
296 * The integral is
|
|
|
297 *
|
|
|
298 * erf(x) = 2/ $(SQRT)($(PI))
|
|
|
299 * $(INTEGRATE 0, x) exp( - $(POWER t, 2)) dt
|
|
|
300 *
|
|
|
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.
|
|
|
303 */
|
|
|
304 real erf(real x)
|
|
|
305 {
|
|
|
306 return std.internal.math.errorfunction.erf(x);
|
|
|
307 }
|
|
|
308
|
|
|
309 /** Complementary error function
|
|
|
310 *
|
|
|
311 * erfc(x) = 1 - erf(x)
|
|
|
312 * = 2/ $(SQRT)($(PI))
|
|
|
313 * $(INTEGRATE x, $(INFIN)) exp( - $(POWER t, 2)) dt
|
|
|
314 *
|
|
|
315 * This function has high relative accuracy for
|
|
|
316 * values of x far from zero. (For values near zero, use erf(x)).
|
|
|
317 */
|
|
|
318 real erfc(real x)
|
|
|
319 {
|
|
|
320 return std.internal.math.errorfunction.erfc(x);
|
|
|
321 }
|
|
|
322
|
|
|
323
|
|
|
324 /** Normal distribution function.
|
|
|
325 *
|
|
|
326 * The normal (or Gaussian, or bell-shaped) distribution is
|
|
|
327 * defined as:
|
|
|
328 *
|
|
|
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))
|
|
|
332 *
|
|
|
333 * To maintain accuracy at values of x near 1.0, use
|
|
|
334 * normalDistribution(x) = 1.0 - normalDistribution(-x).
|
|
|
335 *
|
|
|
336 * References:
|
|
|
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).
|
|
|
340 */
|
|
|
341 real normalDistribution(real x)
|
|
|
342 {
|
|
|
343 return std.internal.math.errorfunction.normalDistributionImpl(x);
|
|
|
344 }
|
|
|
345
|
|
|
346 /** Inverse of Normal distribution function
|
|
|
347 *
|
|
|
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.
|
|
|
351 *
|
|
|
352 * Note: This function is only implemented to 80 bit precision.
|
|
|
353 */
|
|
|
354 real normalDistributionInverse(real p)
|
|
|
355 in {
|
|
|
356 assert(p >= 0.0L && p <= 1.0L, "Domain error");
|
|
|
357 }
|
|
|
358 body
|
|
|
359 {
|
|
|
360 return std.internal.math.errorfunction.normalDistributionInvImpl(p);
|
|
|
361 }
|
