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 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 }