# D-Programming-Language/phobos

### Subversion checkout URL

You can clone with
or
.
Fetching contributors…

Cannot retrieve contributors at this time

358 lines (333 sloc) 9.77 KB
 // Written in the D programming language. /** * Mathematical Special Functions * * The technical term 'Special Functions' includes several families of * transcendental functions, which have important applications in particular * branches of mathematics and physics. * * The gamma and related functions, and the error function are crucial for * mathematical statistics. * The Bessel and related functions arise in problems involving wave propagation * (especially in optics). * Other major categories of special functions include the elliptic integrals * (related to the arc length of an ellipse), and the hypergeometric functions. * * Status: * Many more functions will be added to this module. * The naming convention for the distribution functions (gammaIncomplete, etc) * is not yet finalized and will probably change. * * Macros: * WIKI = Phobos/StdMathSpecial * * TABLE_SV = * Special Values * \$0 * SVH = \$(TR \$(TH \$1) \$(TH \$2)) * SV = \$(TR \$(TD \$1) \$(TD \$2)) * * NAN = \$(RED NAN) * SUP = \$0 * GAMMA = Γ * THETA = θ * INTEGRAL = ∫ * INTEGRATE = \$(BIG ∫\$(SMALL \$1)\$2) * POWER = \$1\$2 * SUB = \$1\$2 * BIGSUM = \$(BIG Σ \$2\$(SMALL \$1)) * CHOOSE = \$(BIG () \$(SMALL \$1)\$(SMALL \$2) \$(BIG )) * PLUSMN = ± * INFIN = ∞ * PLUSMNINF = ±∞ * PI = π * LT = < * GT = > * SQRT = √ * HALF = ½ * * * Copyright: Based on the CEPHES math library, which is * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). * License: \$(WEB www.boost.org/LICENSE_1_0.txt, Boost License 1.0). * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston * Source: \$(PHOBOSSRC std/_mathspecial.d) */ module std.mathspecial; public import std.math; private import std.internal.math.gammafunction; private import std.internal.math.errorfunction; /* *********************************************** * GAMMA AND RELATED FUNCTIONS * * ***********************************************/ pure: nothrow: @safe: @nogc: /** The Gamma function, \$(GAMMA)(x) * * \$(GAMMA)(x) is a generalisation of the factorial function * to real and complex numbers. * Like x!, \$(GAMMA)(x+1) = x * \$(GAMMA)(x). * * Mathematically, if z.re > 0 then * \$(GAMMA)(z) = \$(INTEGRATE 0, \$(INFIN)) \$(POWER t, z-1)\$(POWER e, -t) dt * * \$(TABLE_SV * \$(SVH x, \$(GAMMA)(x) ) * \$(SV \$(NAN), \$(NAN) ) * \$(SV \$(PLUSMN)0.0, \$(PLUSMNINF)) * \$(SV integer > 0, (x-1)! ) * \$(SV integer < 0, \$(NAN) ) * \$(SV +\$(INFIN), +\$(INFIN) ) * \$(SV -\$(INFIN), \$(NAN) ) * ) */ real gamma(real x) { return std.internal.math.gammafunction.gamma(x); } /** Natural logarithm of the gamma function, \$(GAMMA)(x) * * Returns the base e (2.718...) logarithm of the absolute * value of the gamma function of the argument. * * For reals, logGamma is equivalent to log(fabs(gamma(x))). * * \$(TABLE_SV * \$(SVH x, logGamma(x) ) * \$(SV \$(NAN), \$(NAN) ) * \$(SV integer <= 0, +\$(INFIN) ) * \$(SV \$(PLUSMNINF), +\$(INFIN) ) * ) */ real logGamma(real x) { return std.internal.math.gammafunction.logGamma(x); } /** The sign of \$(GAMMA)(x). * * Returns -1 if \$(GAMMA)(x) < 0, +1 if \$(GAMMA)(x) > 0, * \$(NAN) if sign is indeterminate. * * Note that this function can be used in conjunction with logGamma(x) to * evaluate gamma for very large values of x. */ real sgnGamma(real x) { /* Author: Don Clugston. */ if (isNaN(x)) return x; if (x > 0) return 1.0; if (x < -1/real.epsilon) { // Large negatives lose all precision return real.nan; } long n = rndtol(x); if (x == n) { return x == 0 ? copysign(1, x) : real.nan; } return n & 1 ? 1.0 : -1.0; } unittest { assert(sgnGamma(5.0) == 1.0); assert(isNaN(sgnGamma(-3.0))); assert(sgnGamma(-0.1) == -1.0); assert(sgnGamma(-55.1) == 1.0); assert(isNaN(sgnGamma(-real.infinity))); assert(isIdentical(sgnGamma(NaN(0xABC)), NaN(0xABC))); } /** Beta function * * The beta function is defined as * * beta(x, y) = (\$(GAMMA)(x) * \$(GAMMA)(y)) / \$(GAMMA)(x + y) */ real beta(real x, real y) { if ((x+y)> MAXGAMMA) { return exp(logGamma(x) + logGamma(y) - logGamma(x+y)); } else return gamma(x) * gamma(y) / gamma(x+y); } unittest { assert(isIdentical(beta(NaN(0xABC), 4), NaN(0xABC))); assert(isIdentical(beta(2, NaN(0xABC)), NaN(0xABC))); } /** Digamma function * * The digamma function is the logarithmic derivative of the gamma function. * * digamma(x) = d/dx logGamma(x) * * See_Also: \$(LREF logmdigamma), \$(LREF logmdigammaInverse). */ real digamma(real x) { return std.internal.math.gammafunction.digamma(x); } /** Log Minus Digamma function * * logmdigamma(x) = log(x) - digamma(x) * * See_Also: \$(LREF digamma), \$(LREF logmdigammaInverse). */ real logmdigamma(real x) { return std.internal.math.gammafunction.logmdigamma(x); } /** Inverse of the Log Minus Digamma function * * Given y, the function finds x such log(x) - digamma(x) = y. * * See_Also: \$(LREF logmdigamma), \$(LREF digamma). */ real logmdigammaInverse(real x) { return std.internal.math.gammafunction.logmdigammaInverse(x); } /** Incomplete beta integral * * Returns incomplete beta integral of the arguments, evaluated * from zero to x. The regularized incomplete beta function is defined as * * betaIncomplete(a, b, x) = \$(GAMMA)(a + b) / ( \$(GAMMA)(a) \$(GAMMA)(b) ) * * \$(INTEGRATE 0, x) \$(POWER t, a-1)\$(POWER (1-t), b-1) dt * * and is the same as the the cumulative distribution function. * * The domain of definition is 0 <= x <= 1. In this * implementation a and b are restricted to positive values. * The integral from x to 1 may be obtained by the symmetry * relation * * betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x ) * * The integral is evaluated by a continued fraction expansion * or, when b * x is small, by a power series. */ real betaIncomplete(real a, real b, real x ) { return std.internal.math.gammafunction.betaIncomplete(a, b, x); } /** Inverse of incomplete beta integral * * Given y, the function finds x such that * * betaIncomplete(a, b, x) == y * * Newton iterations or interval halving is used. */ real betaIncompleteInverse(real a, real b, real y ) { return std.internal.math.gammafunction.betaIncompleteInv(a, b, y); } /** Incomplete gamma integral and its complement * * These functions are defined by * * gammaIncomplete = ( \$(INTEGRATE 0, x) \$(POWER e, -t) \$(POWER t, a-1) dt )/ \$(GAMMA)(a) * * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) * = (\$(INTEGRATE x, \$(INFIN)) \$(POWER e, -t) \$(POWER t, a-1) dt )/ \$(GAMMA)(a) * * In this implementation both arguments must be positive. * The integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. */ real gammaIncomplete(real a, real x ) in { assert(x >= 0); assert(a > 0); } body { return std.internal.math.gammafunction.gammaIncomplete(a, x); } /** ditto */ real gammaIncompleteCompl(real a, real x ) in { assert(x >= 0); assert(a > 0); } body { return std.internal.math.gammafunction.gammaIncompleteCompl(a, x); } /** Inverse of complemented incomplete gamma integral * * Given a and p, the function finds x such that * * gammaIncompleteCompl( a, x ) = p. */ real gammaIncompleteComplInverse(real a, real p) in { assert(p >= 0 && p <= 1); assert(a > 0); } body { return std.internal.math.gammafunction.gammaIncompleteComplInv(a, p); } /* *********************************************** * ERROR FUNCTIONS & NORMAL DISTRIBUTION * * ***********************************************/ /** Error function * * The integral is * * erf(x) = 2/ \$(SQRT)(\$(PI)) * \$(INTEGRATE 0, x) exp( - \$(POWER t, 2)) dt * * The magnitude of x is limited to about 106.56 for IEEE 80-bit * arithmetic; 1 or -1 is returned outside this range. */ real erf(real x) { return std.internal.math.errorfunction.erf(x); } /** Complementary error function * * erfc(x) = 1 - erf(x) * = 2/ \$(SQRT)(\$(PI)) * \$(INTEGRATE x, \$(INFIN)) exp( - \$(POWER t, 2)) dt * * This function has high relative accuracy for * values of x far from zero. (For values near zero, use erf(x)). */ real erfc(real x) { return std.internal.math.errorfunction.erfc(x); } /** Normal distribution function. * * The normal (or Gaussian, or bell-shaped) distribution is * defined as: * * normalDist(x) = 1/\$(SQRT)(2\$(PI)) \$(INTEGRATE -\$(INFIN), x) exp( - \$(POWER t, 2)/2) dt * = 0.5 + 0.5 * erf(x/sqrt(2)) * = 0.5 * erfc(- x/sqrt(2)) * * To maintain accuracy at values of x near 1.0, use * normalDistribution(x) = 1.0 - normalDistribution(-x). * * References: * \$(LINK http://www.netlib.org/cephes/ldoubdoc.html), * G. Marsaglia, "Evaluating the Normal Distribution", * Journal of Statistical Software 11, (July 2004). */ real normalDistribution(real x) { return std.internal.math.errorfunction.normalDistributionImpl(x); } /** Inverse of Normal distribution function * * Returns the argument, x, for which the area under the * Normal probability density function (integrated from * minus infinity to x) is equal to p. */ real normalDistributionInverse(real p) in { assert(p>=0.0L && p<=1.0L, "Domain error"); } body { return std.internal.math.errorfunction.normalDistributionInvImpl(p); }
Something went wrong with that request. Please try again.