Mercurial > octave-dspies

changeset 8361:cf620941af1a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Set max_recursion_depth and use a subfunction in nchoosek
author Francesco Potortì <pot@gnu.org>
date Sat, 29 Nov 2008 17:57:50 +0100
parents 32d9c0299e99
children 03b414516dd8
files scripts/ChangeLog scripts/specfun/nchoosek.m
diffstat 2 files changed, 24 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Sat Nov 29 00:26:57 2008 +0100
+++ b/scripts/ChangeLog	Sat Nov 29 17:57:50 2008 +0100
@@ -2,6 +2,10 @@
 
 	* plot/__go_draw_axes__.m: Set two point clipping mode to be on.
 
+2008-11-26  Francesco Potortì  <pot@gnu.org>
+
+	* specfun/nchoosek.m: Set max_recursion_depth and use a subfunction.
+
 2008-11-24  Ben Abbott <bpabbott@mac.com>
 
 	* plot/legend.m: Correct ording of legend labels.
--- a/scripts/specfun/nchoosek.m	Sat Nov 29 00:26:57 2008 +0100
+++ b/scripts/specfun/nchoosek.m	Sat Nov 29 17:57:50 2008 +0100
@@ -75,9 +75,26 @@
   elseif (k == n)
      A = v(:).';
   else
-    m = round (exp (sum (log (k:n-1)) - sum (log (2:n-k))));
-    A = [v(1)*ones(m,1), nchoosek(v(2:n),k-1);
-	 nchoosek(v(2:n),k)];
+    oldmax = max_recursion_depth ();
+    unwind_protect
+      max_recursion_depth (n);
+      A = nck (v, k);
+    unwind_protect_cleanup
+      max_recursion_depth (oldmax);
+    end_unwind_protect
   endif
+endfunction
 
+function A = nck (v, k)
+  n = length (v);
+  if (n == 1 || k < 2 || k == n)
+    A = nchoosek (v, k);
+  else
+    m = nchoosek (n-1, k-1);
+    A = [v(1)*ones(m,1), nck(v(2:n),k-1);
+	 nck(v(2:n), k)];
+  endif
 endfunction
+
+%!assert (nchoosek(100,45), bincoeff(100,45))
+%!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])