new changeset for nchoosek recursion problem

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

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ì <[hidden email]>
# 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ì  <pot@gnu.
-2008-11-25  Francesco Potortì  <[hidden email]>
-
- * specfun/nchoosek.m: Conditionally set max_recursion_depth.
+2008-11-26  Francesco Potortì  <[hidden email]>
+
+ * specfun/nchoosek.m: Set max_recursion_depth and use a subfunction.
 
 2008-10-14  Francesco Potortì  <[hidden email]>
 
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: [hidden email]
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/