src/nmath/gamma_cody.c - R - Git at Google

 /* From http://www.netlib.org/specfun/gamma	Fortran translated by f2c,...
  *	------------------------------#####	Martin Maechler, ETH Zurich
  *
  *=========== was part of	ribesl (Bessel I(.))
  *===========			~~~~~~
  */

 // used in bessel_i.c and bessel_j.c, hidden if possible.

 #include "nmath.h"

 double attribute_hidden Rf_gamma_cody(double x)
 {
 /* ----------------------------------------------------------------------

    This routine calculates the GAMMA function for a float argument X.
    Computation is based on an algorithm outlined in reference [1].
    The program uses rational functions that approximate the GAMMA
    function to at least 20 significant decimal digits.	Coefficients
    for the approximation over the interval (1,2) are unpublished.
    Those for the approximation for X >= 12 are from reference [2].
    The accuracy achieved depends on the arithmetic system, the
    compiler, the intrinsic functions, and proper selection of the
    machine-dependent constants.

    *******************************************************************

    Error returns

    The program returns the value XINF for singularities or
    when overflow would occur.	 The computation is believed
    to be free of underflow and overflow.

    Intrinsic functions required are:

    INT, DBLE, EXP, LOG, REAL, SIN


    References:
    [1]  "An Overview of Software Development for Special Functions",
 	W. J. Cody, Lecture Notes in Mathematics, 506,
 	Numerical Analysis Dundee, 1975, G. A. Watson (ed.),
 	Springer Verlag, Berlin, 1976.

    [2]  Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.

    Latest modification: October 12, 1989

    Authors: W. J. Cody and L. Stoltz
    Applied Mathematics Division
    Argonne National Laboratory
    Argonne, IL 60439
    ----------------------------------------------------------------------*/

 /* ----------------------------------------------------------------------
    Mathematical constants
    ----------------------------------------------------------------------*/
     const static double sqrtpi = .9189385332046727417803297; /* == ??? */

 /* *******************************************************************

    Explanation of machine-dependent constants

    beta	- radix for the floating-point representation
    maxexp - the smallest positive power of beta that overflows
    XBIG	- the largest argument for which GAMMA(X) is representable
 	in the machine, i.e., the solution to the equation
 	GAMMA(XBIG) = beta**maxexp
    XINF	- the largest machine representable floating-point number;
 	approximately beta**maxexp
    EPS	- the smallest positive floating-point number such that  1.0+EPS > 1.0
    XMININ - the smallest positive floating-point number such that
 	1/XMININ is machine representable

    Approximate values for some important machines are:

    beta	      maxexp	     XBIG

    CRAY-1		(S.P.)	      2		8191	    966.961
    Cyber 180/855
    under NOS	(S.P.)	      2		1070	    177.803
    IEEE (IBM/XT,
    SUN, etc.)	(S.P.)	      2		 128	    35.040
    IEEE (IBM/XT,
    SUN, etc.)	(D.P.)	      2		1024	    171.624
    IBM 3033	(D.P.)	     16		  63	    57.574
    VAX D-Format	(D.P.)	      2		 127	    34.844
    VAX G-Format	(D.P.)	      2		1023	    171.489

    XINF	 EPS	    XMININ

    CRAY-1		(S.P.)	 5.45E+2465   7.11E-15	  1.84E-2466
    Cyber 180/855
    under NOS	(S.P.)	 1.26E+322    3.55E-15	  3.14E-294
    IEEE (IBM/XT,
    SUN, etc.)	(S.P.)	 3.40E+38     1.19E-7	  1.18E-38
    IEEE (IBM/XT,
    SUN, etc.)	(D.P.)	 1.79D+308    2.22D-16	  2.23D-308
    IBM 3033	(D.P.)	 7.23D+75     2.22D-16	  1.39D-76
    VAX D-Format	(D.P.)	 1.70D+38     1.39D-17	  5.88D-39
    VAX G-Format	(D.P.)	 8.98D+307    1.11D-16	  1.12D-308

    *******************************************************************

    ----------------------------------------------------------------------
    Machine dependent parameters
    ----------------------------------------------------------------------
    */


     const static double xbig = 171.624;
     /* ML_POSINF ==   const double xinf = 1.79e308;*/
     /* DBL_EPSILON = const double eps = 2.22e-16;*/
     /* DBL_MIN ==   const double xminin = 2.23e-308;*/

     /*----------------------------------------------------------------------
       Numerator and denominator coefficients for rational minimax
       approximation over (1,2).
       ----------------------------------------------------------------------*/
     const static double p[8] = {
 	-1.71618513886549492533811,
 	24.7656508055759199108314,-379.804256470945635097577,
 	629.331155312818442661052,866.966202790413211295064,
 	-31451.2729688483675254357,-36144.4134186911729807069,
 	66456.1438202405440627855 };
     const static double q[8] = {
 	-30.8402300119738975254353,
 	315.350626979604161529144,-1015.15636749021914166146,
 	-3107.77167157231109440444,22538.1184209801510330112,
 	4755.84627752788110767815,-134659.959864969306392456,
 	-115132.259675553483497211 };
     /*----------------------------------------------------------------------
       Coefficients for minimax approximation over (12, INF).
       ----------------------------------------------------------------------*/
     const static double c[7] = {
 	-.001910444077728,8.4171387781295e-4,
 	-5.952379913043012e-4,7.93650793500350248e-4,
 	-.002777777777777681622553,.08333333333333333331554247,
 	.0057083835261 };

     /* Local variables */
     int i, n;
     int parity;/*logical*/
     double fact, xden, xnum, y, z, yi, res, sum, ysq;

     parity = (0);
     fact = 1.;
     n = 0;
     y = x;
     if (y <= 0.) {
 	/* -------------------------------------------------------------
 	   Argument is negative
 	   ------------------------------------------------------------- */
 	y = -x;
 	yi = trunc(y);
 	res = y - yi;
 	if (res != 0.) {
 	    if (yi != trunc(yi * .5) * 2.)
 		parity = (1);
 	    fact = -M_PI / sinpi(res);
 	    y += 1.;
 	} else {
 	    return(ML_POSINF);
 	}
     }
     /* -----------------------------------------------------------------
        Argument is positive
        -----------------------------------------------------------------*/
     if (y < DBL_EPSILON) {
 	/* --------------------------------------------------------------
 	   Argument < EPS
 	   -------------------------------------------------------------- */
 	if (y >= DBL_MIN) {
 	    res = 1. / y;
 	} else {
 	    return(ML_POSINF);
 	}
     } else if (y < 12.) {
 	yi = y;
 	if (y < 1.) {
 	    /* ---------------------------------------------------------
 	       EPS < argument < 1
 	       --------------------------------------------------------- */
 	    z = y;
 	    y += 1.;
 	} else {
 	    /* -----------------------------------------------------------
 	       1 <= argument < 12, reduce argument if necessary
 	       ----------------------------------------------------------- */
 	    n = (int) y - 1;
 	    y -= (double) n;
 	    z = y - 1.;
 	}
 	/* ---------------------------------------------------------
 	   Evaluate approximation for 1. < argument < 2.
 	   ---------------------------------------------------------*/
 	xnum = 0.;
 	xden = 1.;
 	for (i = 0; i < 8; ++i) {
 	    xnum = (xnum + p[i]) * z;
 	    xden = xden * z + q[i];
 	}
 	res = xnum / xden + 1.;
 	if (yi < y) {
 	    /* --------------------------------------------------------
 	       Adjust result for case  0. < argument < 1.
 	       -------------------------------------------------------- */
 	    res /= yi;
 	} else if (yi > y) {
 	    /* ----------------------------------------------------------
 	       Adjust result for case  2. < argument < 12.
 	       ---------------------------------------------------------- */
 	    for (i = 0; i < n; ++i) {
 		res *= y;
 		y += 1.;
 	    }
 	}
     } else {
 	/* -------------------------------------------------------------
 	   Evaluate for argument >= 12.,
 	   ------------------------------------------------------------- */
 	if (y <= xbig) {
 	    ysq = y * y;
 	    sum = c[6];
 	    for (i = 0; i < 6; ++i) {
 		sum = sum / ysq + c[i];
 	    }
 	    sum = sum / y - y + sqrtpi;
 	    sum += (y - .5) * log(y);
 	    res = exp(sum);
 	} else {
 	    return(ML_POSINF);
 	}
     }
     /* ----------------------------------------------------------------------
        Final adjustments and return
        ----------------------------------------------------------------------*/
     if (parity)
 	res = -res;
     if (fact != 1.)
 	res = fact / res;
     return res;
 }
	/* From http://www.netlib.org/specfun/gamma Fortran translated by f2c,...
	* ------------------------------##### Martin Maechler, ETH Zurich
	*
	*=========== was part of ribesl (Bessel I(.))
	*=========== ~~~~~~
	*/

	// used in bessel_i.c and bessel_j.c, hidden if possible.

	#include "nmath.h"

	double attribute_hidden Rf_gamma_cody(double x)
	{
	/* ----------------------------------------------------------------------

	This routine calculates the GAMMA function for a float argument X.
	Computation is based on an algorithm outlined in reference [1].
	The program uses rational functions that approximate the GAMMA
	function to at least 20 significant decimal digits. Coefficients
	for the approximation over the interval (1,2) are unpublished.
	Those for the approximation for X >= 12 are from reference [2].
	The accuracy achieved depends on the arithmetic system, the
	compiler, the intrinsic functions, and proper selection of the
	machine-dependent constants.

	*******************************************************************

	Error returns

	The program returns the value XINF for singularities or
	when overflow would occur. The computation is believed
	to be free of underflow and overflow.

	Intrinsic functions required are:

	INT, DBLE, EXP, LOG, REAL, SIN


	References:
	[1] "An Overview of Software Development for Special Functions",
	W. J. Cody, Lecture Notes in Mathematics, 506,
	Numerical Analysis Dundee, 1975, G. A. Watson (ed.),
	Springer Verlag, Berlin, 1976.

	[2] Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.

	Latest modification: October 12, 1989

	Authors: W. J. Cody and L. Stoltz
	Applied Mathematics Division
	Argonne National Laboratory
	Argonne, IL 60439
	----------------------------------------------------------------------*/

	/* ----------------------------------------------------------------------
	Mathematical constants
	----------------------------------------------------------------------*/
	const static double sqrtpi = .9189385332046727417803297; /* == ??? */

	/* *******************************************************************

	Explanation of machine-dependent constants

	beta - radix for the floating-point representation
	maxexp - the smallest positive power of beta that overflows
	XBIG - the largest argument for which GAMMA(X) is representable
	in the machine, i.e., the solution to the equation
	GAMMA(XBIG) = beta**maxexp
	XINF - the largest machine representable floating-point number;
	approximately beta**maxexp
	EPS - the smallest positive floating-point number such that 1.0+EPS > 1.0
	XMININ - the smallest positive floating-point number such that
	1/XMININ is machine representable

	Approximate values for some important machines are:

	beta maxexp XBIG

	CRAY-1 (S.P.) 2 8191 966.961
	Cyber 180/855
	under NOS (S.P.) 2 1070 177.803
	IEEE (IBM/XT,
	SUN, etc.) (S.P.) 2 128 35.040
	IEEE (IBM/XT,
	SUN, etc.) (D.P.) 2 1024 171.624
	IBM 3033 (D.P.) 16 63 57.574
	VAX D-Format (D.P.) 2 127 34.844
	VAX G-Format (D.P.) 2 1023 171.489

	XINF EPS XMININ

	CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
	Cyber 180/855
	under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
	IEEE (IBM/XT,
	SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
	IEEE (IBM/XT,
	SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
	IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
	VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39
	VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308

	*******************************************************************

	----------------------------------------------------------------------
	Machine dependent parameters
	----------------------------------------------------------------------
	*/


	const static double xbig = 171.624;
	/* ML_POSINF == const double xinf = 1.79e308;*/
	/* DBL_EPSILON = const double eps = 2.22e-16;*/
	/* DBL_MIN == const double xminin = 2.23e-308;*/

	/*----------------------------------------------------------------------
	Numerator and denominator coefficients for rational minimax
	approximation over (1,2).
	----------------------------------------------------------------------*/
	const static double p[8] = {
	-1.71618513886549492533811,
	24.7656508055759199108314,-379.804256470945635097577,
	629.331155312818442661052,866.966202790413211295064,
	-31451.2729688483675254357,-36144.4134186911729807069,
	66456.1438202405440627855 };
	const static double q[8] = {
	-30.8402300119738975254353,
	315.350626979604161529144,-1015.15636749021914166146,
	-3107.77167157231109440444,22538.1184209801510330112,
	4755.84627752788110767815,-134659.959864969306392456,
	-115132.259675553483497211 };
	/*----------------------------------------------------------------------
	Coefficients for minimax approximation over (12, INF).
	----------------------------------------------------------------------*/
	const static double c[7] = {
	-.001910444077728,8.4171387781295e-4,
	-5.952379913043012e-4,7.93650793500350248e-4,
	-.002777777777777681622553,.08333333333333333331554247,
	.0057083835261 };

	/* Local variables */
	int i, n;
	int parity;/logical/
	double fact, xden, xnum, y, z, yi, res, sum, ysq;

	parity = (0);
	fact = 1.;
	n = 0;
	y = x;
	if (y <= 0.) {
	/* -------------------------------------------------------------
	Argument is negative
	------------------------------------------------------------- */
	y = -x;
	yi = trunc(y);
	res = y - yi;
	if (res != 0.) {
	if (yi != trunc(yi * .5) * 2.)
	parity = (1);
	fact = -M_PI / sinpi(res);
	y += 1.;
	} else {
	return(ML_POSINF);
	}
	}
	/* -----------------------------------------------------------------
	Argument is positive
	-----------------------------------------------------------------*/
	if (y < DBL_EPSILON) {
	/* --------------------------------------------------------------
	Argument < EPS
	-------------------------------------------------------------- */
	if (y >= DBL_MIN) {
	res = 1. / y;
	} else {
	return(ML_POSINF);
	}
	} else if (y < 12.) {
	yi = y;
	if (y < 1.) {
	/* ---------------------------------------------------------
	EPS < argument < 1
	--------------------------------------------------------- */
	z = y;
	y += 1.;
	} else {
	/* -----------------------------------------------------------
	1 <= argument < 12, reduce argument if necessary
	----------------------------------------------------------- */
	n = (int) y - 1;
	y -= (double) n;
	z = y - 1.;
	}
	/* ---------------------------------------------------------
	Evaluate approximation for 1. < argument < 2.
	---------------------------------------------------------*/
	xnum = 0.;
	xden = 1.;
	for (i = 0; i < 8; ++i) {
	xnum = (xnum + p[i]) * z;
	xden = xden * z + q[i];
	}
	res = xnum / xden + 1.;
	if (yi < y) {
	/* --------------------------------------------------------
	Adjust result for case 0. < argument < 1.
	-------------------------------------------------------- */
	res /= yi;
	} else if (yi > y) {
	/* ----------------------------------------------------------
	Adjust result for case 2. < argument < 12.
	---------------------------------------------------------- */
	for (i = 0; i < n; ++i) {
	res *= y;
	y += 1.;
	}
	}
	} else {
	/* -------------------------------------------------------------
	Evaluate for argument >= 12.,
	------------------------------------------------------------- */
	if (y <= xbig) {
	ysq = y * y;
	sum = c[6];
	for (i = 0; i < 6; ++i) {
	sum = sum / ysq + c[i];
	}
	sum = sum / y - y + sqrtpi;
	sum += (y - .5) * log(y);
	res = exp(sum);
	} else {
	return(ML_POSINF);
	}
	}
	/* ----------------------------------------------------------------------
	Final adjustments and return
	----------------------------------------------------------------------*/
	if (parity)
	res = -res;
	if (fact != 1.)
	res = fact / res;
	return res;
	}