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

annotate scripts/specfun/nchoosek.m @ 10549:95c3e38098bf

Untabify .m scripts

author	Rik <code@nomad.inbox5.com>
date	Fri, 23 Apr 2010 11:28:50 -0700
parents	a7a9eecc07b5
children	23b3ae008f5e

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 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	27 ## @tex
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	28 ## $$
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	29 ## {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	30 ## = {n! \over k! (n-k)!}
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	31 ## $$
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	32 ## @end tex
8517 81d6ab3ac93c Allow documentation tobe built for other formats than tex and info sh@sh-laptop parents: 8506 diff changeset	33 ## @ifnottex
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	34 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	35 ## @example
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	36 ## @group
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	37 ## / \
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	38 ## \| 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	39 ## \| \| = ------------------------- = ---------
451b346d8c2f [project @ 2007-06-25 17:31:46 by jwe] jwe parents: 6575 diff changeset	40 ## \| k \| k! k! (n-k)!
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	41 ## \ /
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	42 ## @end group
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	43 ## @end example
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));
6318 bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	91 else
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	92 p = cell (1, k);
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 8404 diff changeset	93 ## Hack: do the op in the smallest integer class possible to avoid
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 8404 diff changeset	94 ## moving too much data.
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	95 if (n < intmax ("uint8"))
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	96 cl = "uint8";
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	97 elseif (n < intmax ("uint16"))
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	98 cl = "uint16";
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	99 elseif (n < intmax ("uint32"))
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	100 cl = "uint32";
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	101 else
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 8404 diff changeset	102 ## This would exhaust memory anyway.
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	103 cl = "double";
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	104 endif
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	105
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 8404 diff changeset	106 ## Use a generalized Pascal triangle. Traverse backwards to keep
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 8404 diff changeset	107 ## alphabetical order.
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	108 for i = 1:k
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	109 p{i} = zeros (0, i, cl);
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	110 endfor
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	111 s = ones (1, 1, cl);
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	112 p{1} = n*s;
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	113 for j = n-1:-1:1
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	114 for i = k:-1:2
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9210 diff changeset	115 q = p{i-1};
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9210 diff changeset	116 p{i} = [[repmat(s*j, rows (p{i-1}), 1), p{i-1}]; p{i}];
8391 343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	117 endfor
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	118 p{1} = [j;p{1}];
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	119 endfor
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	120 v = v(:);
343f0fbca6eb implement nchoosek without recursion Jaroslav Hajek <highegg@gmail.com> parents: 8361 diff changeset	121 A = v(p{k});
6318 bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	122 endif
8361 cf620941af1a Set max_recursion_depth and use a subfunction in nchoosek Francesco Potortì <pot@gnu.org> parents: 7017 diff changeset	123 endfunction
6318 bae85c1e0e2a [project @ 2007-02-16 08:15:25 by jwe] jwe parents: 6316 diff changeset	124
8404 868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	125 %!warning (nchoosek(100,45));
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	126 %!error (nchoosek(100,45.5));
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	127 %!error (nchoosek(100,145));
868149aac690 nchoosek checks, warning, corner cases and tests Francesco Potortì <pot@gnu.org> parents: 8391 diff changeset	128 %!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	129 %!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 @ 10549:95c3e38098bf