-- Bug description --
There is no nchoosek function in Scilab.
This is a suggestion of implementation, which takes
into account for floating point issues.
I release this code under the Cecill.
The tests are included in the comments.
//
// nchoosek --
// Returns the binomial number (n,k), i.e.
// the number of k-element subsets of an n-element set.
// It is defined by n!/(k! (n-k)!
// and can be computed as
//
// n (n-1) ... (n-k+1)
// -------------------
// k (k-1) ... 1
//
// References
// http://en.wikipedia.org/wiki/Binomial_coefficients
// http://wiki.tcl.tk/1755
//
// Note about floating point accuracy
// factorial(170) ~ 7.e306
// So n=170 is the greatest integer for which
// n! can be computed.
// If the naive formula
//
// b = factorial(n)/factorial(k)/factorial(n-k);
//
// was used, 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 > 1.e308.
// This is why the gammaln function 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 know compute c = b (mod 4) and you
// get c = b = 3.99...99998, intead of getting c = 4.
// This is why, when input arguments are integers,
// the result is rounded to the nearest integer.
// c = nchoosek ( 4 , 1 ) // 4
// c = nchoosek ( 5 , 0 ) // 1
// c = nchoosek ( 5 , 1 ) // 5
// c = nchoosek ( 5 , 2 ) // 10
// c = nchoosek ( 5 , 3 ) // 10
// c = nchoosek ( 5 , 4 ) // 5
// c = nchoosek ( 5 , 5 ) // 1
// c = nchoosek ( 17 , 18 ) // 0
// c = nchoosek ( 17 , -1 ) // 0
// c = nchoosek ( 1.5 , 0.5 ) // 1.5
// c = nchoosek ( 10000 , 134 ) // 2.050083865024873735e307
//
function b = nchoosek ( n , k )
if ( ( k < 0 ) | ( k > n ) ) then
b = 0
else
r = gammaln ( n + 1 ) - gammaln (k + 1) - gammaln (n - k + 1)
b = exp( r );
end
// If the input where integers, returns also an integer.
if ( and(round(n)==n) & and(round(k)==k) ) then
b = round ( b )
end
endfunction
-- Scilab error message --
-- How to reproduce the bug --