# HG changeset patch
# User  Francesco Potortì  <pot@gnu.org>
# Date 1227977870 -3600
# Node ID cf620941af1a8b727ba014cd87f667e97dbfaa69
# Parent  32d9c0299e9921810d52365b1db488f2458cbf07
Set max_recursion_depth and use a subfunction in nchoosek

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -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.
diff --git a/scripts/specfun/nchoosek.m b/scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m
+++ b/scripts/specfun/nchoosek.m
@@ -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])