1######################################################################## 2## 3## Copyright (C) 1995-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {} bincoeff (@var{n}, @var{k}) 28## Return the binomial coefficient of @var{n} and @var{k}. 29## 30## The binomial coefficient is defined as 31## @tex 32## $$ 33## {n \choose k} = {n (n-1) (n-2) \cdots (n-k+1) \over k!} 34## $$ 35## @end tex 36## @ifnottex 37## 38## @example 39## @group 40## / \ 41## | n | n (n-1) (n-2) @dots{} (n-k+1) 42## | | = ------------------------- 43## | k | k! 44## \ / 45## @end group 46## @end example 47## 48## @end ifnottex 49## For example: 50## 51## @example 52## @group 53## bincoeff (5, 2) 54## @result{} 10 55## @end group 56## @end example 57## 58## In most cases, the @code{nchoosek} function is faster for small 59## scalar integer arguments. It also warns about loss of precision for 60## big arguments. 61## 62## @seealso{nchoosek} 63## @end deftypefn 64 65function b = bincoeff (n, k) 66 67 if (nargin != 2) 68 print_usage (); 69 endif 70 71 [retval, n, k] = common_size (n, k); 72 if (retval > 0) 73 error ("bincoeff: N and K must be of common size or scalars"); 74 endif 75 76 if (iscomplex (n) || iscomplex (k)) 77 error ("bincoeff: N and K must not be complex"); 78 endif 79 80 b = zeros (size (n)); 81 82 ok = (k >= 0) & (k == fix (k)) & (! isnan (n)); 83 b(! ok) = NaN; 84 85 n_int = (n == fix (n)); 86 idx = n_int & (n < 0) & ok; 87 b(idx) = (-1) .^ k(idx) .* exp (gammaln (abs (n(idx)) + k(idx)) 88 - gammaln (k(idx) + 1) 89 - gammaln (abs (n(idx)))); 90 91 idx = (n >= k) & ok; 92 b(idx) = exp (gammaln (n(idx) + 1) 93 - gammaln (k(idx) + 1) 94 - gammaln (n(idx) - k(idx) + 1)); 95 96 idx = (! n_int) & (n < k) & ok; 97 b(idx) = (1/pi) * exp (gammaln (n(idx) + 1) 98 - gammaln (k(idx) + 1) 99 + gammaln (k(idx) - n(idx)) 100 + log (sin (pi * (n(idx) - k(idx) + 1)))); 101 102 ## Clean up rounding errors. 103 b(n_int) = round (b(n_int)); 104 105 idx = ! n_int; 106 b(idx) = real (b(idx)); 107 108endfunction 109 110 111%!assert (bincoeff (4, 2), 6) 112%!assert (bincoeff (2, 4), 0) 113%!assert (bincoeff (-4, 2), 10) 114%!assert (bincoeff (5, 2), 10) 115%!assert (bincoeff (50, 6), 15890700) 116%!assert (bincoeff (0.4, 2), -.12, 8*eps) 117 118%!assert (bincoeff ([4 NaN 4], [-1, 2, 2.5]), NaN (1, 3)) 119 120## Test input validation 121%!error bincoeff () 122%!error bincoeff (1, 2, 3) 123%!error bincoeff (ones (3),ones (2)) 124%!error bincoeff (ones (2),ones (3)) 125