| /* |
| * Mathlib : A C Library of Special Functions |
| * Copyright (C) 2000-2018 The R Core Team |
| * Copyright (C) 2002-2018 The R Foundation |
| * Copyright (C) 1998 Ross Ihaka |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| * |
| * SYNOPSIS |
| * |
| * #include <Rmath.h> |
| * double gammafn(double x); |
| * |
| * DESCRIPTION |
| * |
| * This function computes the value of the gamma function. |
| * |
| * NOTES |
| * |
| * This function is a translation into C of a Fortran subroutine |
| * by W. Fullerton of Los Alamos Scientific Laboratory. |
| * (e.g. http://www.netlib.org/slatec/fnlib/gamma.f) |
| * |
| * The accuracy of this routine compares (very) favourably |
| * with those of the Sun Microsystems portable mathematical |
| * library. |
| * |
| * MM specialized the case of n! for n < 50 - for even better precision |
| */ |
| |
| #include "nmath.h" |
| |
| double gammafn(double x) |
| { |
| const static double gamcs[42] = { |
| +.8571195590989331421920062399942e-2, |
| +.4415381324841006757191315771652e-2, |
| +.5685043681599363378632664588789e-1, |
| -.4219835396418560501012500186624e-2, |
| +.1326808181212460220584006796352e-2, |
| -.1893024529798880432523947023886e-3, |
| +.3606925327441245256578082217225e-4, |
| -.6056761904460864218485548290365e-5, |
| +.1055829546302283344731823509093e-5, |
| -.1811967365542384048291855891166e-6, |
| +.3117724964715322277790254593169e-7, |
| -.5354219639019687140874081024347e-8, |
| +.9193275519859588946887786825940e-9, |
| -.1577941280288339761767423273953e-9, |
| +.2707980622934954543266540433089e-10, |
| -.4646818653825730144081661058933e-11, |
| +.7973350192007419656460767175359e-12, |
| -.1368078209830916025799499172309e-12, |
| +.2347319486563800657233471771688e-13, |
| -.4027432614949066932766570534699e-14, |
| +.6910051747372100912138336975257e-15, |
| -.1185584500221992907052387126192e-15, |
| +.2034148542496373955201026051932e-16, |
| -.3490054341717405849274012949108e-17, |
| +.5987993856485305567135051066026e-18, |
| -.1027378057872228074490069778431e-18, |
| +.1762702816060529824942759660748e-19, |
| -.3024320653735306260958772112042e-20, |
| +.5188914660218397839717833550506e-21, |
| -.8902770842456576692449251601066e-22, |
| +.1527474068493342602274596891306e-22, |
| -.2620731256187362900257328332799e-23, |
| +.4496464047830538670331046570666e-24, |
| -.7714712731336877911703901525333e-25, |
| +.1323635453126044036486572714666e-25, |
| -.2270999412942928816702313813333e-26, |
| +.3896418998003991449320816639999e-27, |
| -.6685198115125953327792127999999e-28, |
| +.1146998663140024384347613866666e-28, |
| -.1967938586345134677295103999999e-29, |
| +.3376448816585338090334890666666e-30, |
| -.5793070335782135784625493333333e-31 |
| }; |
| |
| int i, n; |
| double y; |
| double sinpiy, value; |
| |
| #ifdef NOMORE_FOR_THREADS |
| static int ngam = 0; |
| static double xmin = 0, xmax = 0., xsml = 0., dxrel = 0.; |
| |
| /* Initialize machine dependent constants, the first time gamma() is called. |
| FIXME for threads ! */ |
| if (ngam == 0) { |
| ngam = chebyshev_init(gamcs, 42, DBL_EPSILON/20);/*was .1*d1mach(3)*/ |
| gammalims(&xmin, &xmax);/*-> ./gammalims.c */ |
| xsml = exp(fmax2(log(DBL_MIN), -log(DBL_MAX)) + 0.01); |
| /* = exp(.01)*DBL_MIN = 2.247e-308 for IEEE */ |
| dxrel = sqrt(DBL_EPSILON);/*was sqrt(d1mach(4)) */ |
| } |
| #else |
| /* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 : |
| * (xmin, xmax) are non-trivial, see ./gammalims.c |
| * xsml = exp(.01)*DBL_MIN |
| * dxrel = sqrt(DBL_EPSILON) = 2 ^ -26 |
| */ |
| # define ngam 22 |
| # define xmin -170.5674972726612 |
| # define xmax 171.61447887182298 |
| # define xsml 2.2474362225598545e-308 |
| # define dxrel 1.490116119384765696e-8 |
| #endif |
| |
| if(ISNAN(x)) return x; |
| |
| /* If the argument is exactly zero or a negative integer |
| * then return NaN. */ |
| if (x == 0 || (x < 0 && x == round(x))) { |
| ML_ERROR(ME_DOMAIN, "gammafn"); |
| return ML_NAN; |
| } |
| |
| y = fabs(x); |
| |
| if (y <= 10) { |
| |
| /* Compute gamma(x) for -10 <= x <= 10 |
| * Reduce the interval and find gamma(1 + y) for 0 <= y < 1 |
| * first of all. */ |
| |
| n = (int) x; |
| if(x < 0) --n; |
| y = x - n;/* n = floor(x) ==> y in [ 0, 1 ) */ |
| --n; |
| value = chebyshev_eval(y * 2 - 1, gamcs, ngam) + .9375; |
| if (n == 0) |
| return value;/* x = 1.dddd = 1+y */ |
| |
| if (n < 0) { |
| /* compute gamma(x) for -10 <= x < 1 */ |
| |
| /* exact 0 or "-n" checked already above */ |
| |
| /* The answer is less than half precision */ |
| /* because x too near a negative integer. */ |
| if (x < -0.5 && fabs(x - (int)(x - 0.5) / x) < dxrel) { |
| ML_ERROR(ME_PRECISION, "gammafn"); |
| } |
| |
| /* The argument is so close to 0 that the result would overflow. */ |
| if (y < xsml) { |
| ML_ERROR(ME_RANGE, "gammafn"); |
| if(x > 0) return ML_POSINF; |
| else return ML_NEGINF; |
| } |
| |
| n = -n; |
| |
| for (i = 0; i < n; i++) { |
| value /= (x + i); |
| } |
| return value; |
| } |
| else { |
| /* gamma(x) for 2 <= x <= 10 */ |
| |
| for (i = 1; i <= n; i++) { |
| value *= (y + i); |
| } |
| return value; |
| } |
| } |
| else { |
| /* gamma(x) for y = |x| > 10. */ |
| |
| if (x > xmax) { /* Overflow */ |
| // No warning: +Inf is the best answer |
| return ML_POSINF; |
| } |
| |
| if (x < xmin) { /* Underflow */ |
| // No warning: 0 is the best answer |
| return 0.; |
| } |
| |
| if(y <= 50 && y == (int)y) { /* compute (n - 1)! */ |
| value = 1.; |
| for (i = 2; i < y; i++) value *= i; |
| } |
| else { /* normal case */ |
| value = exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI + |
| ((2*y == (int)2*y)? stirlerr(y) : lgammacor(y))); |
| } |
| if (x > 0) |
| return value; |
| |
| if (fabs((x - (int)(x - 0.5))/x) < dxrel){ |
| |
| /* The answer is less than half precision because */ |
| /* the argument is too near a negative integer. */ |
| |
| ML_ERROR(ME_PRECISION, "gammafn"); |
| } |
| |
| sinpiy = sinpi(y); |
| if (sinpiy == 0) { /* Negative integer arg - overflow */ |
| ML_ERROR(ME_RANGE, "gammafn"); |
| return ML_POSINF; |
| } |
| |
| return -M_PI / (y * sinpiy * value); |
| } |
| } |