Mercurial > hg > octave-nkf

comparison scripts/specfun/nchoosek.m @ 10800:23b3ae008f5e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
optimize nchoosek
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 19 Jul 2010 13:35:42 +0200
parents 95c3e38098bf
children 693e22af08ae
comparison
equal deleted inserted replaced
10799:177f1ad7c7c1 10800:23b3ae008f5e
86 A = v(:); 86 A = v(:);
87 elseif (k == n) 87 elseif (k == n)
88 A = v(:).'; 88 A = v(:).';
89 elseif (k > n) 89 elseif (k > n)
90 A = zeros (0, k, class (v)); 90 A = zeros (0, k, class (v));
91 else 91 elseif (k == 2)
92 p = cell (1, k); 92 ## Can do it without transpose.
93 ## Hack: do the op in the smallest integer class possible to avoid 93 x = repelems (v(1:n-1), [1:n-1; n-1:-1:1]).';
94 ## moving too much data. 94 y = cat (1, cellslices (v(:), 2:n, n*ones (1, n-1)){:});
95 if (n < intmax ("uint8")) 95 A = [x, y];
96 cl = "uint8"; 96 elseif (k < n)
97 elseif (n < intmax ("uint16")) 97 v = v(:).';
98 cl = "uint16"; 98 A = v(k:n);
99 elseif (n < intmax ("uint32")) 99 l = 1:n-k+1;
100 cl = "uint32"; 100 for j = 2:k
101 else 101 c = columns (A);
102 ## This would exhaust memory anyway. 102 cA = cellslices (A, l, c*ones (1, n-k+1), 2);
103 cl = "double"; 103 l = c-l+1;
104 endif 104 b = repelems (v(k-j+1:n-j+1), [1:n-k+1; l]);
105 105 A = [b; cA{:}];
106 ## Use a generalized Pascal triangle. Traverse backwards to keep 106 l = cumsum (l);
107 ## alphabetical order. 107 l = [1, 1 + l(1:n-k)];
108 for i = 1:k
109 p{i} = zeros (0, i, cl);
110 endfor 108 endfor
111 s = ones (1, 1, cl); 109 clear cA b;
112 p{1} = n*s; 110 A = A.';
113 for j = n-1:-1:1
114 for i = k:-1:2
115 q = p{i-1};
116 p{i} = [[repmat(s*j, rows (p{i-1}), 1), p{i-1}]; p{i}];
117 endfor
118 p{1} = [j;p{1}];
119 endfor
120 v = v(:);
121 A = v(p{k});
122 endif 111 endif
123 endfunction 112 endfunction
124 113
125 %!warning (nchoosek(100,45)); 114 %!warning (nchoosek(100,45));
126 %!error (nchoosek(100,45.5)); 115 %!error (nchoosek(100,45.5));