changelog shortlog tags changeset files revisions annotate raw

scripts/miscellaneous/bincoeff.m

changeset 10289: 4b124317dc38
parent:1bf0ce0930be
author: John W. Eaton <jwe@octave.org>
date: Tue Feb 09 20:58:55 2010 -0500 (47 minutes ago)
permissions: -rw-r--r--
description: base_properties::set_children: account for hidden children 
     1## Copyright (C) 1995, 1996, 1999, 2000, 2002, 2004, 2005, 2006, 2007,

     2##               2008, 2009 Kurt Hornik

     3##

     4## This file is part of Octave.

     5##

     6## Octave is free software; you can redistribute it and/or modify it

     7## under the terms of the GNU General Public License as published by

     8## the Free Software Foundation; either version 3 of the License, or (at

     9## your option) any later version.

    10##

    11## Octave is distributed in the hope that it will be useful, but

    12## WITHOUT ANY WARRANTY; without even the implied warranty of

    13## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

    14## General Public License for more details.

    15##

    16## You should have received a copy of the GNU General Public License

    17## along with Octave; see the file COPYING.  If not, see

    18## <http://www.gnu.org/licenses/>.

    19

    20## -*- texinfo -*-

    21## @deftypefn {Mapping Function} {} bincoeff (@var{n}, @var{k})

    22## Return the binomial coefficient of @var{n} and @var{k}, defined as

    23## @tex

    24## $$

    25##  {n \choose k} = {n (n-1) (n-2) \cdots (n-k+1) \over k!}

    26## $$

    27## @end tex

    28## @ifnottex

    29##

    30## @example

    31## @group

    32##  /   \

    33##  | n |    n (n-1) (n-2) @dots{} (n-k+1)

    34##  |   |  = -------------------------

    35##  | k |               k!

    36##  \   /

    37## @end group

    38## @end example

    39## @end ifnottex

    40##

    41## For example,

    42##

    43## @example

    44## @group

    45## bincoeff (5, 2)

    46##      @result{} 10

    47## @end group

    48## @end example

    49##

    50## In most cases, the @code{nchoosek} function is faster for small

    51## scalar integer arguments.  It also warns about loss of precision for

    52## big arguments.

    53##

    54## @seealso{nchoosek}

    55## @end deftypefn

    56

    57## Author: KH <Kurt.Hornik@wu-wien.ac.at>

    58## Created: 8 October 1994

    59## Adapted-By: jwe

    60

    61function b = bincoeff (n, k)

    62

    63  if (nargin != 2)

    64    print_usage ();

    65  endif

    66

    67  [retval, n, k] = common_size (n, k);

    68  if (retval > 0)

    69    error ("bincoeff: n and k must be of common size or scalars");

    70  endif

    71

    72  sz = size (n);

    73  b   = zeros (sz);

    74

    75  ind = (! (k >= 0) | (k != real (round (k))) | isnan (n));

    76  b(ind) = NaN;

    77  

    78  ind = (k == 0);

    79  b(ind) = 1;

    80

    81  ind = ((k > 0) & ((n == real (round (n))) & (n < 0)));

    82  b(ind) = (-1) .^ k(ind) .* exp (gammaln (abs (n(ind)) + k(ind))

    83				  - gammaln (k(ind) + 1)

    84				  - gammaln (abs (n(ind))));

    85

    86  ind = ((k > 0) & (n >= k));

    87  b(ind) = exp (gammaln (n(ind) + 1)

    88		- gammaln (k(ind) + 1)

    89		- gammaln (n(ind) - k(ind) + 1));

    90

    91  ind = ((k > 0) & ((n != real (round (n))) & (n < k)));

    92  b(ind) = (1/pi) * exp (gammaln (n(ind) + 1)

    93			 - gammaln (k(ind) + 1)

    94			 + gammaln (k(ind) - n(ind))

    95			 + log (sin (pi * (n(ind) - k(ind) + 1))));

    96

    97  ## Clean up rounding errors.

    98  ind = (n == round (n));

    99  b(ind) = round (b(ind));

   100

   101  ind = (n != round (n));

   102  b(ind) = real (b(ind));

   103

   104endfunction

   105

   106%!assert(bincoeff(4,2), 6)

   107%!assert(bincoeff(2,4), 0)

   108%!assert(bincoeff(0.4,2), -.12, 8*eps)

   109

   110%!assert(bincoeff (5, 2) == 10 && bincoeff (50, 6) == 15890700);

   111

   112%!error bincoeff ();

   113

   114%!error bincoeff (1, 2, 3);