rmathd.choose source code

1 module rmathd.choose;
2 
3 public import rmathd.common;
4 public import rmathd.beta;
5 
6 
7 /*
8 ** dmd choose.d common.d beta.d && ./choose
9 */
10 
11 
12 T lfastchoose(T)(T n, T k)
13 {
14     return -log(n + 1.) - lbeta!T(n - k + 1., k + 1.);
15 }
16 
17 
18 /* mathematically the same:
19    less stable typically, but useful if n-k+1 < 0 : */
20 static T lfastchoose2(T)(T n, T k, int *s_choose)
21 {
22     T r;
23     r = lgammafn_sign(n - k + 1., s_choose);
24     return lgammafn(n + 1.) - lgammafn(k + 1.) - r;
25 }
26 
27 
28 template ODD(T){
29 	T ODD(T k){
30 		return (k != 2 * floor(k/2.0));
31 	}
32 }
33 
34 
35 T lchoose(T: double)(T n, T k)
36 {
37     T k0 = k;
38     k = nearbyint(k);
39     
40     /* NaNs propagated correctly */
41     if(isNaN(n) || isNaN(k))
42     	return n + k;
43     
44     if (fabs(k - k0) > 1e-7){
45 	    //MATHLIB_WARNING2(_("'k' (%.2f) must be integer, rounded to %.0f"), k0, k);
46 	    doNothing();
47     }
48     
49     if (k < 2) {
50 	    if (k <	 0)
51 	    	return -T.infinity;
52 	    if (k == 0)
53 	    	return 0.;
54 	    /* else: k == 1 */
55 	    return log(fabs(n));
56     }
57     /* else: k >= 2 */
58     if (n < 0) {
59 	    return lchoose!T(-n + k - 1, k);
60     } else if (!R_nonint(n)) {
61 	    n = nearbyint(n);
62 	    if(n < k)
63 	    	return -T.infinity;
64 	    /* k <= n :*/
65 	    if(n - k < 2)
66 	    	return lchoose!T(n, n - k); /* <- Symmetry */
67 	    /* else: n >= k+2 */
68 	    return lfastchoose!T(n, k);
69     }
70     /* else non-integer n >= 0 : */
71     if (n < k - 1) {
72 	    int s;
73 	    return lfastchoose2!T(n, k, &s);
74     }
75     return lfastchoose!T(n, k);
76 }
77 
78 
79 enum k_small_max = 30;
80 /* 30 is somewhat arbitrary: it is on the *safe* side:
81  * both speed and precision are clearly improved for k < 30.
82 */
83 
84 
85 T choose(T: double)(T n, T k)
86 {
87     T r, k0 = k;
88     k = nearbyint(k);
89     
90     /* NaNs propagated correctly */
91     if(isNaN(n) || isNaN(k))
92     	return n + k;
93     
94     if (fabs(k - k0) > 1e-7){
95 	    //MATHLIB_WARNING2(_("'k' (%.2f) must be integer, rounded to %.0f"), k0, k);
96 	    doNothing();
97     }
98 
99     if (k < k_small_max) {
100 	    int j;
101 	    if(n - k < k && n >= 0 && !R_nonint(n))
102 	    	k = n - k; /* <- Symmetry */
103 	    if (k <	 0)
104 	    	return 0.;
105 	    if (k == 0)
106 	    	return 1.;
107 	    /* else: k >= 1 */
108 	    r = n;
109 	    for(j = 2; j <= k; j++)
110 	        r *= (n - j + 1)/j;
111 	    return !R_nonint(n) ? nearbyint(r) : r;
112 	    /* might have got rounding errors */
113     }
114     /* else: k >= k_small_max */
115     if (n < 0) {
116 	    r = choose!T(-n + k - 1, k);
117 	    if (ODD!T(k)) r = -r;
118 	    return r;
119     }
120     else if (!R_nonint(n)) {
121 	    n = nearbyint(n);
122 	    if(n < k)
123 	    	return 0.;
124 	    if(n - k < k_small_max)
125 	    	return choose!T(n, n - k); /* <- Symmetry */
126 	    return nearbyint(exp(lfastchoose!T(n, k)));
127     }
128     /* else non-integer n >= 0 : */
129     if (n < k - 1) {
130 	    int s_choose;
131 	    r = lfastchoose2(n, k, /* -> */ &s_choose);
132 	    return s_choose * exp(r);
133     }
134     return exp(lfastchoose!T(n, k));
135 }
136 
137 
138 void test_choose()
139 {
140 	import std.stdio: writeln;
141 	writeln("choose(7., 4.): ", choose(7., 4.));
142 	writeln("lchoose(7., 4.): ", lchoose(7., 4.));
143 }
144