The R Project SVN R

Rev

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

Rev 11736 Rev 24736
Line 38... Line 38...
38 
#include "nmath.h"
 38 
#include "nmath.h"
39 
#include "dpq.h"
 39 
#include "dpq.h"
40 
 
 40 
 
41 
double dgamma(double x, double shape, double scale, int give_log)
 41 
double dgamma(double x, double shape, double scale, int give_log)
42 
{
42 
{
43 
#ifndef D_non_pois
 - 
 
44 
    double pr;
 43 
    double pr;
45 
#endif
 - 
 
46 
#ifdef IEEE_754
 44 
#ifdef IEEE_754
47 
    if (ISNAN(x) || ISNAN(shape) || ISNAN(scale))
 45 
    if (ISNAN(x) || ISNAN(shape) || ISNAN(scale))
48 
        return x + shape + scale;
 46 
        return x + shape + scale;
49 
#endif
 47 
#endif
50 
    if (shape <= 0 || scale <= 0) ML_ERR_return_NAN;
 48 
    if (shape <= 0 || scale <= 0) ML_ERR_return_NAN;
Line 55... Line 53...
55 
	if (shape > 1) return R_D__0;
 53 
	if (shape > 1) return R_D__0;
56 
	/* else */
 54 
	/* else */
57 
	return give_log ? -log(scale) : 1 / scale;
 55 
	return give_log ? -log(scale) : 1 / scale;
58 
    }
 56 
    }
59
57
60 
#ifdef D_non_pois
 - 
 
61 
 
 - 
 
62 
    x /= scale;
 - 
 
63 
    return give_log ?
 - 
 
64 
	   ((shape - 1) * log(x) - lgammafn(shape) - x) - log(scale) :
 - 
 
65 
	exp((shape - 1) * log(x) - lgammafn(shape) - x) / scale;
 - 
 
66 
 
 - 
 
67 
#else /* new dpois() based code */
 - 
 
68 
 
 - 
 
69 
    if (shape < 1) {
58 
    if (shape < 1) {
70 
	pr = dpois_raw(shape, x/scale, give_log);
 59 
	pr = dpois_raw(shape, x/scale, give_log);
71 
	return give_log ?  pr + log(shape/x) : pr*shape/x;
 60 
	return give_log ?  pr + log(shape/x) : pr*shape/x;
72 
    }
 61 
    }
73 
    /* else  shape >= 1 */
 62 
    /* else  shape >= 1 */
74 
    pr = dpois_raw(shape-1, x/scale, give_log);
 63 
    pr = dpois_raw(shape-1, x/scale, give_log);
75 
    return give_log ? pr - log(scale) : pr/scale;
 64 
    return give_log ? pr - log(scale) : pr/scale;
76 
#endif
 - 
 
77 
}
 65 
}