The R Project SVN R

Rev

Rev 77685 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 77685 Rev 85878
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) 2000-2018 The R Core Team
 3 
 *  Copyright (C) 2000-2024 The R Core Team
4 
 *  Copyright (C) 2002-2018 The R Foundation
 4 
 *  Copyright (C) 2002-2024 The R Foundation
5 
 *  Copyright (C) 1998 Ross Ihaka
 5 
 *  Copyright (C) 1998 Ross Ihaka
6 
 *
 6 
 *
7 
 *  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
8 
 *  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
9 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
Line 87... Line 87...
87 
	-.1967938586345134677295103999999e-29,
 87 
	-.1967938586345134677295103999999e-29,
88 
	+.3376448816585338090334890666666e-30,
 88 
	+.3376448816585338090334890666666e-30,
89 
	-.5793070335782135784625493333333e-31
 89 
	-.5793070335782135784625493333333e-31
90 
    };
 90 
    };
91 
 
 91 
 
92 
    int i, n;
 - 
 
93 
    double y;
 - 
 
94 
    double sinpiy, value;
 - 
 
95 
 
 - 
 
96 
#ifdef NOMORE_FOR_THREADS
 92 
#ifdef NOMORE_FOR_THREADS
97 
    static int ngam = 0;
 93 
    static int ngam = 0;
98 
    static double xmin = 0, xmax = 0., xsml = 0., dxrel = 0.;
 94 
    static double xmin = 0, xmax = 0., xsml = 0., dxrel = 0.;
99 
 
 95 
 
100 
    /* Initialize machine dependent constants, the first time gamma() is called.
 96 
    /* Initialize machine dependent constants, the first time gamma() is called.
Line 126... Line 122...
126 
    if (x == 0 || (x < 0 && x == round(x))) {
 122 
    if (x == 0 || (x < 0 && x == round(x))) {
127 
	ML_WARNING(ME_DOMAIN, "gammafn");
 123 
	ML_WARNING(ME_DOMAIN, "gammafn");
128 
	return ML_NAN;
 124 
	return ML_NAN;
129 
    }
 125 
    }
130 
 
 126 
 
- 
 
 127 
    int i;
131 
    y = fabs(x);
 128 
    double y = fabs(x), value;
132 
 
 129 
 
133 
    if (y <= 10) {
 130 
    if (y <= 10) {
134 
 
 131 
 
135 
	/* Compute gamma(x) for -10 <= x <= 10
 132 
	/* Compute gamma(x) for -10 <= x <= 10
136 
	 * Reduce the interval and find gamma(1 + y) for 0 <= y < 1
 133 
	 * Reduce the interval and find gamma(1 + y) for 0 <= y < 1
137 
	 * first of all. */
 134 
	 * first of all. */
138 
 
 135 
 
139 
	n = (int) x;
 136 
	int n = (int) x;
140 
	if(x < 0) --n;
 137 
	if(x < 0) --n;
141 
	y = x - n;/* n = floor(x)  ==>	y in [ 0, 1 ) */
 138 
	y = x - n;/* n = floor(x)  ==>	y in [ 0, 1 ) */
142 
	--n;
 139 
	--n;
143 
	value = chebyshev_eval(y * 2 - 1, gamcs, ngam) + .9375;
 140 
	value = chebyshev_eval(y * 2 - 1, gamcs, ngam) + .9375;
144 
	if (n == 0)
 141 
	if (n == 0)
Line 195... Line 192...
195 
	    value = 1.;
 192 
	    value = 1.;
196 
	    for (i = 2; i < y; i++) value *= i;
 193 
	    for (i = 2; i < y; i++) value *= i;
197 
	}
 194 
	}
198 
	else { /* normal case */
 195 
	else { /* normal case */
199 
	    value = exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
 196 
	    value = exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
200 
			((2*y == (int)2*y)? stirlerr(y) : lgammacor(y)));
 197 
			((2*y == (int)2*y) ? stirlerr(y) : lgammacor(y)));
201 
	}
 198 
	}
- 
 
 199 
 
202 
	if (x > 0)
 200 
	if (x > 0)
203 
	    return value;
 201 
	    return value;
- 
 
 202 
	// else:  x < 0, not an integer :
204 
 
 203 
 
205 
	if (fabs((x - (int)(x - 0.5))/x) < dxrel){
 204 
	if (fabs((x - (int)(x - 0.5))/x) < dxrel) {
206 
 
 - 
 
207 
	    /* The answer is less than half precision because */
 205 
	    /* The answer is less than half precision because */
208 
	    /* the argument is too near a negative integer. */
 206 
	    /* the argument is too near a negative integer. */
209 
 
 207 
 
210 
	    ML_WARNING(ME_PRECISION, "gammafn");
 208 
	    ML_WARNING(ME_PRECISION, "gammafn");
211 
	}
 209 
	}
212 
 
 210 
 
213 
	sinpiy = sinpi(y);
 211 
	double sinpiy = sinpi(y);
214 
	if (sinpiy == 0) {		/* Negative integer arg - overflow */
 212 
	if (sinpiy == 0) {		/* Negative integer arg - overflow */
215 
	    ML_WARNING(ME_RANGE, "gammafn");
 213 
	    ML_WARNING(ME_RANGE, "gammafn");
216 
	    return ML_POSINF;
 214 
	    return ML_POSINF;
217 
	}
 215 
	}
218 
 
 216