The R Project SVN
R
(root)/branches/R-BIGSEXP/src/nmath/gamma_cody.c – Rev 68553
Rev 68552 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
/* 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 GAMMAfunction to at least 20 significant decimal digits. Coefficientsfor 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, thecompiler, the intrinsic functions, and proper selection of themachine-dependent constants.*******************************************************************Error returnsThe program returns the value XINF for singularities orwhen overflow would occur. The computation is believedto be free of underflow and overflow.Intrinsic functions required are:INT, DBLE, EXP, LOG, REAL, SINReferences:[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, 1989Authors: W. J. Cody and L. StoltzApplied Mathematics DivisionArgonne National LaboratoryArgonne, IL 60439----------------------------------------------------------------------*//* ----------------------------------------------------------------------Mathematical constants----------------------------------------------------------------------*/const static double sqrtpi = .9189385332046727417803297; /* == ??? *//* *******************************************************************Explanation of machine-dependent constantsbeta - radix for the floating-point representationmaxexp - the smallest positive power of beta that overflowsXBIG - the largest argument for which GAMMA(X) is representablein the machine, i.e., the solution to the equationGAMMA(XBIG) = beta**maxexpXINF - the largest machine representable floating-point number;approximately beta**maxexpEPS - the smallest positive floating-point number such that 1.0+EPS > 1.0XMININ - the smallest positive floating-point number such that1/XMININ is machine representableApproximate values for some important machines are:beta maxexp XBIGCRAY-1 (S.P.) 2 8191 966.961Cyber 180/855under NOS (S.P.) 2 1070 177.803IEEE (IBM/XT,SUN, etc.) (S.P.) 2 128 35.040IEEE (IBM/XT,SUN, etc.) (D.P.) 2 1024 171.624IBM 3033 (D.P.) 16 63 57.574VAX D-Format (D.P.) 2 127 34.844VAX G-Format (D.P.) 2 1023 171.489XINF EPS XMININCRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466Cyber 180/855under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294IEEE (IBM/XT,SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38IEEE (IBM/XT,SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39VAX 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 minimaxapproximation 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;}