nchoosek() is available on any pocket calculator for undergraduate students, but not natively in Scilab. Please fix that.
#### Reported by **Michael BAUDIN** ``` -- 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 -- ```
issue