libphobos/src/std/mathspecial.d

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