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 */ 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 /* 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 */ 154real 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 */ 176real 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 */ 187real 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 */ 198real 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 */ 223real 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 */ 236real 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 */ 255real gammaIncomplete(real a, real x ) 256in { 257 assert(x >= 0); 258 assert(a > 0); 259} 260body { 261 return std.internal.math.gammafunction.gammaIncomplete(a, x); 262} 263 264/** ditto */ 265real gammaIncompleteCompl(real a, real x ) 266in { 267 assert(x >= 0); 268 assert(a > 0); 269} 270body { 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 */ 280real gammaIncompleteComplInverse(real a, real p) 281in { 282 assert(p >= 0 && p <= 1); 283 assert(a > 0); 284} 285body { 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 */ 304real 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 */ 318real 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 */ 341real 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 */ 354real normalDistributionInverse(real p) 355in { 356 assert(p >= 0.0L && p <= 1.0L, "Domain error"); 357} 358body 359{ 360 return std.internal.math.errorfunction.normalDistributionInvImpl(p); 361} 362