The R Project SVN R

Rev

Rev 42302 | Rev 64358 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 42302 Rev 42303
Line -... Line 1...
- 
 
 1 
/* From http://www.netlib.org/specfun/gamma	Fortran translated by f2c,...
- 
 
 2 
 *	------------------------------#####	Martin Maechler, ETH Zurich
- 
 
 3 
 *
- 
 
 4 
 *=========== was part of	ribesl (Bessel I(.))
- 
 
 5 
 *===========			~~~~~~
- 
 
 6 
 */
- 
 
 7 
#include "nmath.h"
- 
 
 8 
 
- 
 
 9 
double attribute_hidden gamma_cody(double x)
- 
 
 10 
{
- 
 
 11 
/* ----------------------------------------------------------------------
- 
 
 12 
 
- 
 
 13 
   This routine calculates the GAMMA function for a float argument X.
- 
 
 14 
   Computation is based on an algorithm outlined in reference [1].
- 
 
 15 
   The program uses rational functions that approximate the GAMMA
- 
 
 16 
   function to at least 20 significant decimal digits.	Coefficients
- 
 
 17 
   for the approximation over the interval (1,2) are unpublished.
- 
 
 18 
   Those for the approximation for X >= 12 are from reference [2].
- 
 
 19 
   The accuracy achieved depends on the arithmetic system, the
- 
 
 20 
   compiler, the intrinsic functions, and proper selection of the
- 
 
 21 
   machine-dependent constants.
- 
 
 22 
 
- 
 
 23 
   *******************************************************************
- 
 
 24 
 
- 
 
 25 
   Error returns
- 
 
 26 
 
- 
 
 27 
   The program returns the value XINF for singularities or
- 
 
 28 
   when overflow would occur.	 The computation is believed
- 
 
 29 
   to be free of underflow and overflow.
- 
 
 30 
 
- 
 
 31 
   Intrinsic functions required are:
- 
 
 32 
 
- 
 
 33 
   INT, DBLE, EXP, LOG, REAL, SIN
- 
 
 34 
 
- 
 
 35 
 
- 
 
 36 
   References:
- 
 
 37 
   [1]  "An Overview of Software Development for Special Functions",
- 
 
 38 
	W. J. Cody, Lecture Notes in Mathematics, 506,
- 
 
 39 
	Numerical Analysis Dundee, 1975, G. A. Watson (ed.),
- 
 
 40 
	Springer Verlag, Berlin, 1976.
- 
 
 41 
 
- 
 
 42 
   [2]  Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
- 
 
 43 
 
- 
 
 44 
   Latest modification: October 12, 1989
- 
 
 45 
 
- 
 
 46 
   Authors: W. J. Cody and L. Stoltz
- 
 
 47 
   Applied Mathematics Division
- 
 
 48 
   Argonne National Laboratory
- 
 
 49 
   Argonne, IL 60439
- 
 
 50 
   ----------------------------------------------------------------------*/
- 
 
 51 
 
- 
 
 52 
/* ----------------------------------------------------------------------
- 
 
 53 
   Mathematical constants
- 
 
 54 
   ----------------------------------------------------------------------*/
- 
 
 55 
    const static double sqrtpi = .9189385332046727417803297; /* == ??? */
- 
 
 56 
 
- 
 
 57 
/* *******************************************************************
- 
 
 58 
 
- 
 
 59 
   Explanation of machine-dependent constants
- 
 
 60 
 
- 
 
 61 
   beta	- radix for the floating-point representation
- 
 
 62 
   maxexp - the smallest positive power of beta that overflows
- 
 
 63 
   XBIG	- the largest argument for which GAMMA(X) is representable
- 
 
 64 
	in the machine, i.e., the solution to the equation
- 
 
 65 
	GAMMA(XBIG) = beta**maxexp
- 
 
 66 
   XINF	- the largest machine representable floating-point number;
- 
 
 67 
	approximately beta**maxexp
- 
 
 68 
   EPS	- the smallest positive floating-point number such that  1.0+EPS > 1.0
- 
 
 69 
   XMININ - the smallest positive floating-point number such that
- 
 
 70 
	1/XMININ is machine representable
- 
 
 71 
 
- 
 
 72 
   Approximate values for some important machines are:
- 
 
 73 
 
- 
 
 74 
   beta	      maxexp	     XBIG
- 
 
 75 
 
- 
 
 76 
   CRAY-1		(S.P.)	      2		8191	    966.961
- 
 
 77 
   Cyber 180/855
- 
 
 78 
   under NOS	(S.P.)	      2		1070	    177.803
- 
 
 79 
   IEEE (IBM/XT,
- 
 
 80 
   SUN, etc.)	(S.P.)	      2		 128	    35.040
- 
 
 81 
   IEEE (IBM/XT,
- 
 
 82 
   SUN, etc.)	(D.P.)	      2		1024	    171.624
- 
 
 83 
   IBM 3033	(D.P.)	     16		  63	    57.574
- 
 
 84 
   VAX D-Format	(D.P.)	      2		 127	    34.844
- 
 
 85 
   VAX G-Format	(D.P.)	      2		1023	    171.489
- 
 
 86 
 
- 
 
 87 
   XINF	 EPS	    XMININ
- 
 
 88 
 
- 
 
 89 
   CRAY-1		(S.P.)	 5.45E+2465   7.11E-15	  1.84E-2466
- 
 
 90 
   Cyber 180/855
- 
 
 91 
   under NOS	(S.P.)	 1.26E+322    3.55E-15	  3.14E-294
- 
 
 92 
   IEEE (IBM/XT,
- 
 
 93 
   SUN, etc.)	(S.P.)	 3.40E+38     1.19E-7	  1.18E-38
- 
 
 94 
   IEEE (IBM/XT,
- 
 
 95 
   SUN, etc.)	(D.P.)	 1.79D+308    2.22D-16	  2.23D-308
- 
 
 96 
   IBM 3033	(D.P.)	 7.23D+75     2.22D-16	  1.39D-76
- 
 
 97 
   VAX D-Format	(D.P.)	 1.70D+38     1.39D-17	  5.88D-39
- 
 
 98 
   VAX G-Format	(D.P.)	 8.98D+307    1.11D-16	  1.12D-308
- 
 
 99 
 
- 
 
 100 
   *******************************************************************
- 
 
 101 
 
- 
 
 102 
   ----------------------------------------------------------------------
- 
 
 103 
   Machine dependent parameters
- 
 
 104 
   ----------------------------------------------------------------------
- 
 
 105 
   */
- 
 
 106 
 
- 
 
 107 
 
- 
 
 108 
    const static double xbig = 171.624;
- 
 
 109 
    /* ML_POSINF ==   const double xinf = 1.79e308;*/
- 
 
 110 
    /* DBL_EPSILON = const double eps = 2.22e-16;*/
- 
 
 111 
    /* DBL_MIN ==   const double xminin = 2.23e-308;*/
- 
 
 112 
 
- 
 
 113 
    /*----------------------------------------------------------------------
- 
 
 114 
      Numerator and denominator coefficients for rational minimax
- 
 
 115 
      approximation over (1,2).
- 
 
 116 
      ----------------------------------------------------------------------*/
- 
 
 117 
    const static double p[8] = {
- 
 
 118 
	-1.71618513886549492533811,
- 
 
 119 
	24.7656508055759199108314,-379.804256470945635097577,
- 
 
 120 
	629.331155312818442661052,866.966202790413211295064,
- 
 
 121 
	-31451.2729688483675254357,-36144.4134186911729807069,
- 
 
 122 
	66456.1438202405440627855 };
- 
 
 123 
    const static double q[8] = {
- 
 
 124 
	-30.8402300119738975254353,
- 
 
 125 
	315.350626979604161529144,-1015.15636749021914166146,
- 
 
 126 
	-3107.77167157231109440444,22538.1184209801510330112,
- 
 
 127 
	4755.84627752788110767815,-134659.959864969306392456,
- 
 
 128 
	-115132.259675553483497211 };
- 
 
 129 
    /*----------------------------------------------------------------------
- 
 
 130 
      Coefficients for minimax approximation over (12, INF).
- 
 
 131 
      ----------------------------------------------------------------------*/
- 
 
 132 
    const static double c[7] = {
- 
 
 133 
	-.001910444077728,8.4171387781295e-4,
- 
 
 134 
	-5.952379913043012e-4,7.93650793500350248e-4,
- 
 
 135 
	-.002777777777777681622553,.08333333333333333331554247,
- 
 
 136 
	.0057083835261 };
- 
 
 137 
 
- 
 
 138 
    /* Local variables */
- 
 
 139 
    long i, n;
- 
 
 140 
    long int parity;/*logical*/
- 
 
 141 
    double fact, xden, xnum, y, z, yi, res, sum, ysq;
- 
 
 142 
 
- 
 
 143 
    parity = (0);
- 
 
 144 
    fact = 1.;
- 
 
 145 
    n = 0;
- 
 
 146 
    y = x;
- 
 
 147 
    if (y <= 0.) {
- 
 
 148 
	/* -------------------------------------------------------------
- 
 
 149 
	   Argument is negative
- 
 
 150 
	   ------------------------------------------------------------- */
- 
 
 151 
	y = -x;
- 
 
 152 
	yi = ftrunc(y);
- 
 
 153 
	res = y - yi;
- 
 
 154 
	if (res != 0.) {
- 
 
 155 
	    if (yi != ftrunc(yi * .5) * 2.)
- 
 
 156 
		parity = (1);
- 
 
 157 
	    fact = -M_PI / sin(M_PI * res);
- 
 
 158 
	    y += 1.;
- 
 
 159 
	} else {
- 
 
 160 
	    return(ML_POSINF);
- 
 
 161 
	}
- 
 
 162 
    }
- 
 
 163 
    /* -----------------------------------------------------------------
- 
 
 164 
       Argument is positive
- 
 
 165 
       -----------------------------------------------------------------*/
- 
 
 166 
    if (y < DBL_EPSILON) {
- 
 
 167 
	/* --------------------------------------------------------------
- 
 
 168 
	   Argument < EPS
- 
 
 169 
	   -------------------------------------------------------------- */
- 
 
 170 
	if (y >= DBL_MIN) {
- 
 
 171 
	    res = 1. / y;
- 
 
 172 
	} else {
- 
 
 173 
	    return(ML_POSINF);
- 
 
 174 
	}
- 
 
 175 
    } else if (y < 12.) {
- 
 
 176 
	yi = y;
- 
 
 177 
	if (y < 1.) {
- 
 
 178 
	    /* ---------------------------------------------------------
- 
 
 179 
	       EPS < argument < 1
- 
 
 180 
	       --------------------------------------------------------- */
- 
 
 181 
	    z = y;
- 
 
 182 
	    y += 1.;
- 
 
 183 
	} else {
- 
 
 184 
	    /* -----------------------------------------------------------
- 
 
 185 
	       1 <= argument < 12, reduce argument if necessary
- 
 
 186 
	       ----------------------------------------------------------- */
- 
 
 187 
	    n = (long) y - 1;
- 
 
 188 
	    y -= (double) n;
- 
 
 189 
	    z = y - 1.;
- 
 
 190 
	}
- 
 
 191 
	/* ---------------------------------------------------------
- 
 
 192 
	   Evaluate approximation for 1. < argument < 2.
- 
 
 193 
	   ---------------------------------------------------------*/
- 
 
 194 
	xnum = 0.;
- 
 
 195 
	xden = 1.;
- 
 
 196 
	for (i = 0; i < 8; ++i) {
- 
 
 197 
	    xnum = (xnum + p[i]) * z;
- 
 
 198 
	    xden = xden * z + q[i];
- 
 
 199 
	}
- 
 
 200 
	res = xnum / xden + 1.;
- 
 
 201 
	if (yi < y) {
- 
 
 202 
	    /* --------------------------------------------------------
- 
 
 203 
	       Adjust result for case  0. < argument < 1.
- 
 
 204 
	       -------------------------------------------------------- */
- 
 
 205 
	    res /= yi;
- 
 
 206 
	} else if (yi > y) {
- 
 
 207 
	    /* ----------------------------------------------------------
- 
 
 208 
	       Adjust result for case  2. < argument < 12.
- 
 
 209 
	       ---------------------------------------------------------- */
- 
 
 210 
	    for (i = 0; i < n; ++i) {
- 
 
 211 
		res *= y;
- 
 
 212 
		y += 1.;
- 
 
 213 
	    }
- 
 
 214 
	}
- 
 
 215 
    } else {
- 
 
 216 
	/* -------------------------------------------------------------
- 
 
 217 
	   Evaluate for argument >= 12.,
- 
 
 218 
	   ------------------------------------------------------------- */
- 
 
 219 
	if (y <= xbig) {
- 
 
 220 
	    ysq = y * y;
- 
 
 221 
	    sum = c[6];
- 
 
 222 
	    for (i = 0; i < 6; ++i) {
- 
 
 223 
		sum = sum / ysq + c[i];
- 
 
 224 
	    }
- 
 
 225 
	    sum = sum / y - y + sqrtpi;
- 
 
 226 
	    sum += (y - .5) * log(y);
- 
 
 227 
	    res = exp(sum);
- 
 
 228 
	} else {
- 
 
 229 
	    return(ML_POSINF);
- 
 
 230 
	}
- 
 
 231 
    }
- 
 
 232 
    /* ----------------------------------------------------------------------
- 
 
 233 
       Final adjustments and return
- 
 
 234 
       ----------------------------------------------------------------------*/
- 
 
 235 
    if (parity)
- 
 
 236 
	res = -res;
- 
 
 237 
    if (fact != 1.)
- 
 
 238 
	res = fact / res;
- 
 
 239 
    return res;
- 
 
 240 
}
- 
 
 241