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