The R Project SVN R

Rev

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

Rev 19500 Rev 24344
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 1998 Ross Ihaka
 3 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000-2001 The R Development Core Team
 4 
 *  Copyright (C) 2000-2001 The R Development Core Team
- 
 
 5 
 *  Copyright (C) 2002-2003 The R Foundation
5 
 *
 6 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 7 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 8 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 10 
 *  (at your option) any later version.
Line 32... Line 33...
32 
 *    by W. Fullerton of Los Alamos Scientific Laboratory.
 33 
 *    by W. Fullerton of Los Alamos Scientific Laboratory.
33 
 *
 34 
 *
34 
 *    The accuracy of this routine compares (very) favourably
 35 
 *    The accuracy of this routine compares (very) favourably
35 
 *    with those of the Sun Microsystems portable mathematical
 36 
 *    with those of the Sun Microsystems portable mathematical
36 
 *    library.
 37 
 *    library.
- 
 
 38 
 *
- 
 
 39 
 *    MM specialized the case of  n!  for n < 50 - for even better precision
37 
 */
 40 
 */
38 
 
 41 
 
39 
#include "nmath.h"
 42 
#include "nmath.h"
40 
 
 43 
 
41 
double gammafn(double x)
 44 
double gammafn(double x)
Line 102... Line 105...
102 
	/*   = exp(.01)*DBL_MIN = 2.247e-308 for IEEE */
 105 
	/*   = exp(.01)*DBL_MIN = 2.247e-308 for IEEE */
103 
	dxrel = sqrt(1/DBL_EPSILON);/*was (1/d1mach(4)) */
 106 
	dxrel = sqrt(1/DBL_EPSILON);/*was (1/d1mach(4)) */
104 
    }
 107 
    }
105 
#else
 108 
#else
106 
/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
 109 
/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
107 
 * (xmin, xmax) are non-trivial, see ./gammalims.c
110 
 * (xmin, xmax) are non-trivial, see ./gammalims.c
108 
 * xsml = exp(.01)*DBL_MIN
111 
 * xsml = exp(.01)*DBL_MIN
109 
 * dxrel = sqrt(1/DBL_EPSILON) = 2 ^ 26
 112 
 * dxrel = sqrt(1/DBL_EPSILON) = 2 ^ 26
110 
*/
 113 
*/
111 
# define ngam 22
 114 
# define ngam 22
112 
# define xmin -170.5674972726612
 115 
# define xmin -170.5674972726612
113 
# define xmax  171.61447887182298
 116 
# define xmax  171.61447887182298
Line 183... Line 186...
183 
	if (x < xmin) {			/* Underflow */
 186 
	if (x < xmin) {			/* Underflow */
184 
	    ML_ERROR(ME_UNDERFLOW);
 187 
	    ML_ERROR(ME_UNDERFLOW);
185 
	    return ML_UNDERFLOW;
 188 
	    return ML_UNDERFLOW;
186 
	}
 189 
	}
187 
 
 190 
 
- 
 
 191 
	if(y <= 50 && y == (int)y) { /* compute (n - 1)! */
- 
 
 192 
	    value = 1.;
- 
 
 193 
	    for (i = 2; i < y; i++) value *= i;
- 
 
 194 
	}
- 
 
 195 
	else { /* normal case */
188 
	value = exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI + lgammacor(y));
 196 
	    value = exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
- 
 
 197 
			((2*y == (int)2*y)? stirlerr(y) : lgammacor(y)));
189 
 
 198 
	}
190 
	if (x > 0)
 199 
	if (x > 0)
191 
	    return value;
 200 
	    return value;
192 
 
 201 
 
193 
	if (fabs((x - (int)(x - 0.5))/x) < dxrel){
 202 
	if (fabs((x - (int)(x - 0.5))/x) < dxrel){
194 
 
 203