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