specfun_nchoosek — Returns the binomial number (n,k).
a matrix of floating point integers, must be positive
a matrix of floating point integers, must be positive
a matrix of floating point integers
b = specfun_nchoosek ( n , k )
Computes the number of k-element subsets of an n-element set. It is mathematically defined by
which is equivalent to ( n*(n-1)*...*(n-k+1) ) / ( k*(k-1)*...*1 ). The implementation, though, uses a robust floating point implementation, based on the gammaln and exp functions.
If n < 0, k < 0 or k > n, then an error is generated.
Note about floating point accuracy.
We have factorial(170) ~ 7.e306. So n=170 is the greatest integer for which n! can be computed. If the naive formula b = n!/(k! (n-k)! ) was used directly, the maximum value n for which the binomial can be computed is n = 170. But binomial(171,1) = 1, so there is no reason to prevent the computation of 1 because an intermediate result is greater than 1.e308. This is why the gammaln function, combined with exp, is used instead.
Note about rounding for integer inputs.
If n = 4 and k = 1, the gammaln and exp functions are accurate only to 1ulp. This leads to the result that b = 3.99...99998. It is the result of the fact that we use the elementary functions exp and gammaln. This is very close to 4, but is not equal to 4. Assume that you compute c = b (mod 4) and you get c = b = 3.99...99998, intead of getting c = 0. This is why, when input arguments are integers, the result is rounded to the nearest integer with the round function.
c = specfun_nchoosek ( 4 , 1 ) // 4 c = specfun_nchoosek ( 5 , 0 ) // 1 c = specfun_nchoosek ( 5 , 1 ) // 5 c = specfun_nchoosek ( 5 , 2 ) // 10 c = specfun_nchoosek ( 5 , 3 ) // 10 c = specfun_nchoosek ( 5 , 4 ) // 5 c = specfun_nchoosek ( 5 , 5 ) // 1 // The following test shows that the implementation is not naive c = specfun_nchoosek ( 10000 , 134 ) // exact = 2.050083865033972676e307 relerror = abs(c-exact)/exact // On a Windows system, we got ~ 1.e-11, which shows that the implementation // was, in this case, accurate up to 11 digits (instead of 15-17). // The following test shows that the implementation is vectorized. specfun_nchoosek (10,0:10) // [1,10,45,120,210,252,210,120,45,10,1] // Generates an error c = specfun_nchoosek ( 17 , 18 ) c = specfun_nchoosek ( 17 , -1 ) c = specfun_nchoosek ( 1.5 , 0.5 )