new changeset for nchoosek recursion problem

[Francesco Potortì]

Here is a new changeset for nchoosek.  This is against the previous
commit that I did, which worked but was less elegant and, in principle
at least, less efficient.  This changeset also adds tests to nchoosek.


# HG changeset patch
# User Francesco Potortì <address@hidden>
# Date 1227710601 -3600
# Node ID ed030a109bb8bfbeb4def31ac84eebbb975ffa30
# Parent  aed592b1eda936b05eb977135ec8c9bae10a11c0
Set max_recursion_depth and use a subfunction in nchoosek

diff -r aed592b1eda9 -r ed030a109bb8 scripts/ChangeLog
--- a/scripts/ChangeLog Tue Nov 25 15:21:49 2008 +0100
+++ b/scripts/ChangeLog Wed Nov 26 15:43:21 2008 +0100
@@ -1,6 +1,6 @@ 2008-11-25  Francesco Potortì  <address@hidden
-2008-11-25  Francesco Potortì  <address@hidden>
-
-       * specfun/nchoosek.m: Conditionally set max_recursion_depth.
+2008-11-26  Francesco Potortì  <address@hidden>
+
+       * specfun/nchoosek.m: Set max_recursion_depth and use a subfunction.
 
 2008-10-14  Francesco Potortì  <address@hidden>
 
diff -r aed592b1eda9 -r ed030a109bb8 scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m        Tue Nov 25 15:21:49 2008 +0100
+++ b/scripts/specfun/nchoosek.m        Wed Nov 26 15:43:21 2008 +0100
@@ -77,19 +77,24 @@ function A = nchoosek (v, k)
   else
     oldmax = max_recursion_depth ();
     unwind_protect
-      if (n > oldmax)
-       max_recursion_depth (n);
-      else
-       oldmax = false;
-      endif
-      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)];
+      max_recursion_depth (n);
+      A = nck (v, k);
     unwind_protect_cleanup
-      if (oldmax)
-       max_recursion_depth (oldmax);
-      endif
+      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])

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/