The R Project SVN R

Rev

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

Rev 19500 Rev 27260
Line 31... Line 31...
31 
 */
 31 
 */
32 
#include "nmath.h"
 32 
#include "nmath.h"
33 
 
 33 
 
34 
double lfastchoose(double n, double k)
 34 
double lfastchoose(double n, double k)
35 
{
 35 
{
- 
 
 36 
	return -log(n + 1.) - lbeta(n - k + 1., k + 1.);
- 
 
 37 
	/* the same (but less stable):
36 
	return lgammafn(n + 1.0) - lgammafn(k + 1.0) - lgammafn(n - k + 1.0);
 38 
	 * == lgammafn(n + 1.0) - lgammafn(k + 1.0) - lgammafn(n - k + 1.0); */
37 
}
 39 
}
38 
 
 40 
 
39 
double fastchoose(double n, double k)
 41 
double fastchoose(double n, double k)
40 
{
 42 
{
41 
	return exp(lfastchoose(n, k));
 43 
	return exp(lfastchoose(n, k));
Line 47... Line 49...
47 
	k = floor(k + 0.5);
 49 
	k = floor(k + 0.5);
48 
#ifdef IEEE_754
 50 
#ifdef IEEE_754
49 
	/* NaNs propagated correctly */
 51 
	/* NaNs propagated correctly */
50 
	if(ISNAN(n) || ISNAN(k)) return n + k;
 52 
	if(ISNAN(n) || ISNAN(k)) return n + k;
51 
#endif
 53 
#endif
52 
	if (k < 0 || n < k) ML_ERR_return_NAN;
 54 
	if (n < 0) ML_ERR_return_NAN;
- 
 
 55 
	if (k < 0 || n < k) return ML_NEGINF;
53 
 
 56 
 
54 
	return lfastchoose(n, k);
 57 
	return lfastchoose(n, k);
55 
}
 58 
}
56 
 
 59 
 
57 
double choose(double n, double k)
 60 
double choose(double n, double k)
Line 60... Line 63...
60 
	k = floor(k + 0.5);
 63 
	k = floor(k + 0.5);
61 
#ifdef IEEE_754
 64 
#ifdef IEEE_754
62 
	/* NaNs propagated correctly */
 65 
	/* NaNs propagated correctly */
63 
	if(ISNAN(n) || ISNAN(k)) return n + k;
 66 
	if(ISNAN(n) || ISNAN(k)) return n + k;
64 
#endif
 67 
#endif
- 
 
 68 
	if (n < 0) ML_ERR_return_NAN;/* could be defined instead as
- 
 
 69 
					(-1)^k (-n + k - 1 \\ k) {k>=0}*/
65 
	if (k < 0 || n < k) ML_ERR_return_NAN;
 70 
	if (k < 0 || n < k) return 0.;
66 
 
 71 
 
67 
	return floor(exp(lfastchoose(n, k)) + 0.5);
 72 
	return floor(exp(lfastchoose(n, k)) + 0.5);
68 
}
 73 
}