octave_packages/m/miscellaneous/bincoeff.m

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