octave-lyh: scripts/specfun/nchoosek.m annotate

annotate scripts/specfun/nchoosek.m @ 11188:4cb1522e4d0f

Use function handle as input to cellfun, rather than quoted function name or anonymous function wrapper.

author	Rik <octave@nomad.inbox5.com>
date	Wed, 03 Nov 2010 17:20:56 -0700
parents	693e22af08ae
children	fd0a3ac60b0e

rev	line source
8920 eb63fbe60fab update copyright notices John W. Eaton <jwe@octave.org> parents: 8517 diff changeset	1 ## Copyright (C) 2001, 2006, 2007, 2009 Rolf Fabian and Paul Kienzle
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	2 ## Copyright (C) 2008 Jaroslav Hajek
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	3 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	4 ## This file is part of Octave.
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	5 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	6 ## Octave is free software; you can redistribute it and/or modify it
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	7 ## under the terms of the GNU General Public License as published by
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6754 diff changeset	8 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6754 diff changeset	9 ## your option) any later version.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	10 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	11 ## Octave is distributed in the hope that it will be useful, but
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	14 ## General Public License for more details.
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	15 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	16 ## You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6754 diff changeset	17 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6754 diff changeset	18 ## <http://www.gnu.org/licenses/>.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	19
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	20 ## -- texinfo --
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	21 ## @deftypefn {Function File} {@var{c} =} nchoosek (@var{n}, @var{k})
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	22 ##
6575 9ef0450c873d [project @ 2007-04-25 14:07:20 by jwe] jwe parents: 6377 diff changeset	23 ## Compute the binomial coefficient or all combinations of @var{n}.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	24 ## If @var{n} is a scalar then, calculate the binomial coefficient
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	25 ## of @var{n} and @var{k}, defined as
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	26 ## @tex
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	27 ## $$
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	28 ## {n \choose k} = {n (n-1) (n-2) \cdots (n-k+1) \over k!}
6754 451b346d8c2f [project @ 2007-06-25 17:31:46 by jwe] jwe parents: 6575 diff changeset	29 ## = {n! \over k! (n-k)!}
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	30 ## $$
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	31 ## @end tex
8517 81d6ab3ac93c Allow documentation tobe built for other formats than tex and info sh@sh-laptop parents: 8506 diff changeset	32 ## @ifnottex
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	33 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	34 ## @example
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	35 ## @group
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	36 ## / \
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	37 ## \| n \| n (n-1) (n-2) @dots{} (n-k+1) n!
6754 451b346d8c2f [project @ 2007-06-25 17:31:46 by jwe] jwe parents: 6575 diff changeset	38 ## \| \| = ------------------------- = ---------
451b346d8c2f [project @ 2007-06-25 17:31:46 by jwe] jwe parents: 6575 diff changeset	39 ## \| k \| k! k! (n-k)!
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	40 ## \ /
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	41 ## @end group
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	42 ## @end example
10821 693e22af08ae Grammarcheck documentation of m-files Rik <octave@nomad.inbox5.com> parents: 10800 diff changeset	43 ##
8517 81d6ab3ac93c Allow documentation tobe built for other formats than tex and info sh@sh-laptop parents: 8506 diff changeset	44 ## @end ifnottex
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	45 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	46 ## If @var{n} is a vector generate all combinations of the elements
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	47 ## of @var{n}, taken @var{k} at a time, one row per combination. The
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	48 ## resulting @var{c} has size @code{[nchoosek (length (@var{n}),
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	49 ## @var{k}), @var{k}]}.
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	50 ##
9209 923c7cb7f13f Simplify TeXinfo files by eliminating redundant @iftex followed by @tex construction. Rik <rdrider0-list@yahoo.com> parents: 9051 diff changeset	51 ## @code{nchoosek} works only for non-negative integer arguments; use
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	52 ## @code{bincoeff} for non-integer scalar arguments and for using vector
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	53 ## arguments to compute many coefficients at once.
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	54 ##
6316 a3a2580435c2 [project @ 2007-02-16 07:23:49 by jwe] jwe parents: 5827 diff changeset	55 ## @seealso{bincoeff}
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	56 ## @end deftypefn
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	57
6316 a3a2580435c2 [project @ 2007-02-16 07:23:49 by jwe] jwe parents: 5827 diff changeset	58 ## Author: Rolf Fabian <fabian@tu-cottbus.de>
a3a2580435c2 [project @ 2007-02-16 07:23:49 by jwe] jwe parents: 5827 diff changeset	59 ## Author: Paul Kienzle <pkienzle@users.sf.net>
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	60 ## Author: Jaroslav Hajek
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	61
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	62 function A = nchoosek (v, k)
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	63
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	64 if (nargin != 2
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	65 \|\| !isnumeric(k) \|\| !isnumeric(v)
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	66 \|\| !isscalar(k) \|\| (!isscalar(v) && !isvector(v)))
6316 a3a2580435c2 [project @ 2007-02-16 07:23:49 by jwe] jwe parents: 5827 diff changeset	67 print_usage ();
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	68 endif
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	69 if ((isscalar(v) && v < k) \|\| k < 0
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	70 \|\| k != round(k) \|\| any (v < 0 \|\| v != round(v)))
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	71 error ("nchoosek: args are nonnegative integers with V not less than K");
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	72 endif
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	73
6318 bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	74 n = length (v);
bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	75
bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	76 if (n == 1)
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 8404 diff changeset	77 ## Improve precision at next step.
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 8404 diff changeset	78 k = min (k, v-k);
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	79 A = round (prod ((v-k+1:v)./(1:k)));
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	80 if (A2k*eps >= 0.5)
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	81 warning ("nchoosek", "nchoosek: possible loss of precision");
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	82 endif
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	83 elseif (k == 0)
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	84 A = [];
6318 bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	85 elseif (k == 1)
bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	86 A = v(:);
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	87 elseif (k == n)
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	88 A = v(:).';
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	89 elseif (k > n)
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	90 A = zeros (0, k, class (v));
10800 23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	91 elseif (k == 2)
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	92 ## Can do it without transpose.
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	93 x = repelems (v(1:n-1), [1:n-1; n-1:-1:1]).';
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	94 y = cat (1, cellslices (v(:), 2:n, n*ones (1, n-1)){:});
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	95 A = [x, y];
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	96 elseif (k < n)
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	97 v = v(:).';
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	98 A = v(k:n);
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	99 l = 1:n-k+1;
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	100 for j = 2:k
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	101 c = columns (A);
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	102 cA = cellslices (A, l, c*ones (1, n-k+1), 2);
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	103 l = c-l+1;
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	104 b = repelems (v(k-j+1:n-j+1), [1:n-k+1; l]);
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	105 A = [b; cA{:}];
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	106 l = cumsum (l);
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	107 l = [1, 1 + l(1:n-k)];
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	108 endfor
10800 23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	109 clear cA b;
23b3ae008f5e optimize nchoosek Jaroslav Hajek <highegg@gmail.com> parents: 10549 diff changeset	110 A = A.';
6318 bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	111 endif
8361 cf620941af1a Set max_recursion_depth and use a subfunction in nchoosek Francesco Potortì <pot@gnu.org> parents: 7017 diff changeset	112 endfunction
6318 bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	113
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	114 %!warning (nchoosek(100,45));
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	115 %!error (nchoosek(100,45.5));
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	116 %!error (nchoosek(100,145));
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	117 %!assert (nchoosek(80,10), bincoeff(80,10))
8361 cf620941af1a Set max_recursion_depth and use a subfunction in nchoosek Francesco Potortì <pot@gnu.org> parents: 7017 diff changeset	118 %!assert (nchoosek(1:5,3),[1:3;1,2,4;1,2,5;1,3,4;1,3,5;1,4,5;2:4;2,3,5;2,4,5;3:5])

Mercurial > hg > octave-lyh

annotate scripts/specfun/nchoosek.m @ 11188:4cb1522e4d0f