# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1279539342 -7200
# Node ID 23b3ae008f5e7e50906c6822da1dd659c7559b1e
# Parent  177f1ad7c7c14285303c039ac4fa6c3a2837a886
optimize nchoosek

diff -r 177f1ad7c7c1 -r 23b3ae008f5e scripts/ChangeLog
--- a/scripts/ChangeLog	Mon Jul 19 11:59:47 2010 +0200
+++ b/scripts/ChangeLog	Mon Jul 19 13:35:42 2010 +0200
@@ -1,3 +1,7 @@
+2010-07-19  Jaroslav Hajek  <highegg@gmail.com>
+
+	* specfun/nchoosek.m: Optimize.
+
 2010-07-18  Rik <octave@nomad.inbox5.com>
 
 	* scripts/@ftp/ftp.m, scripts/@ftp/mget.m, scripts/audio/lin2mu.m,
diff -r 177f1ad7c7c1 -r 23b3ae008f5e scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m	Mon Jul 19 11:59:47 2010 +0200
+++ b/scripts/specfun/nchoosek.m	Mon Jul 19 13:35:42 2010 +0200
@@ -88,37 +88,26 @@
     A = v(:).';
   elseif (k > n)
     A = zeros (0, k, class (v));
-  else
-    p = cell (1, k);
-    ## Hack: do the op in the smallest integer class possible to avoid
-    ## moving too much data.
-    if (n < intmax ("uint8"))
-      cl = "uint8";
-    elseif (n < intmax ("uint16"))
-      cl = "uint16";
-    elseif (n < intmax ("uint32"))
-      cl = "uint32";
-    else
-      ## This would exhaust memory anyway.
-      cl = "double";
-    endif
-     
-    ## Use a generalized Pascal triangle. Traverse backwards to keep
-    ## alphabetical order.
-    for i = 1:k
-      p{i} = zeros (0, i, cl);
+  elseif (k == 2)
+    ## Can do it without transpose.
+    x = repelems (v(1:n-1), [1:n-1; n-1:-1:1]).';
+    y = cat (1, cellslices (v(:), 2:n, n*ones (1, n-1)){:});
+    A = [x, y];
+  elseif (k < n)
+    v = v(:).';
+    A = v(k:n);
+    l = 1:n-k+1;
+    for j = 2:k
+      c = columns (A);
+      cA = cellslices (A, l, c*ones (1, n-k+1), 2);
+      l = c-l+1;
+      b = repelems (v(k-j+1:n-j+1), [1:n-k+1; l]);
+      A = [b; cA{:}];
+      l = cumsum (l);
+      l = [1, 1 + l(1:n-k)];
     endfor
-    s = ones (1, 1, cl);
-    p{1} = n*s;
-    for j = n-1:-1:1
-      for i = k:-1:2
-        q = p{i-1};
-        p{i} = [[repmat(s*j, rows (p{i-1}), 1), p{i-1}]; p{i}];
-      endfor
-      p{1} = [j;p{1}];
-    endfor
-    v = v(:);
-    A = v(p{k});
+    clear cA b;
+    A = A.';
   endif
 endfunction