1 /** Fixed Length Sorting via Sorting Networks.
2 
3 	See_Also: http://forum.dlang.org/post/ne5m62$1gu5$1@digitalmars.com
4 	See_Also: http://cpansearch.perl.org/src/JGAMBLE/Algorithm-Networksort-1.30/lib/Algorithm/Networksort.pm
5 	See_Also: http://www.angelfire.com/blog/ronz/Articles/999SortingNetworksReferen.html
6 
7 	License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0).
8 
9 	TODO: see some sizes are not supported, we should not have holes.
10 	Use http://www.angelfire.com/blog/ronz/Articles/999SortingNetworksReferen.html
11 
12 	TODO: Sometimes the sort routine gets too bulky. Suggestion: also define
13 	networks for `medianOfUpTo` and `medianExactly`, then use them in a
14 	quicksort manner - first compute the median to segregate values in
15 	below/over the median, then make two calls to `sortExactly!(n / 2)`. That
16 	way you get to sort n values with median of n values (smaller and simpler)
17 	and two sorts of n / 2 values.
18 
19 	TODO: Stability of equal elements: Need template parameter `equalityStability`? Scalar builtin values are always stable.
20 
21 	TODO: There should be a notion of at what point the networks become too bulky
22 	to be fast - 6-16 may be the limit.
23 
24 	TODO: There is a nice peephole optimization you could make. Consider:
25 
26 	r.conditionalSwap!("a < b", less, 0, 1, 1, 2)(r);
27 
28 	which needs to do
29 
30 	if (r[1] < r[0]) r.swapAt(0, 1);
31 	if (r[2] < r[1]) r.swapAt(1, 2);
32 
33 	For types with no elaborate copy/assignment, it's more efficient to use a
34 	"hole"-based approach - assign the first element to a temporary and then
35 	consider it a hole that you fill, then leaving another hole:
36 
37 	if (r[1] < r[0])
38 	if (r[2] < r[0]) r.swapAt(0, 1, 2);
39 	else r.swapAt(0, 1);
40 	else
41 	if (r[2] < r[1]) r.swapAt(1, 2);
42 
43 	with swapAt with three argument having this definition:
44 
45 	auto t = r[a]; r[a] = r[b]; r[b] = r[c]; r[c] = t;
46 
47 	i.e. create a temporary (which creates a "hole" in the array) then fill it
48 	leaving another hole etc., until the last hole is filled with the temporary.
49 */
50 module nxt.algorithm.sortn;
51 
52 import std.meta : allSatisfy;
53 import std.traits : isIntegral;
54 import std.range.primitives : isRandomAccessRange;
55 
56 /** Conditionally pairwise sort elements of `Range` `r` at `indexes` using
57 	comparison predicate `less`.
58 
59 	TODO: Perhaps defines as
60 
61 	template conditionalSwap(indexes...) {
62 	   void conditionalSwap(less = "a < b", Range)(Range r) {
63 	   }
64 	}
65 
66 	instead.
67  */
68 private void conditionalSwap(alias less = "a < b", Range, indexes...)(Range r)
69 if (isRandomAccessRange!Range &&
70 	allSatisfy!(isIntegral, typeof(indexes)) &&
71 	indexes.length &&
72 	(indexes.length & 1) == 0) // even number of indexes
73 {
74 	import std.algorithm.mutation : swapAt;
75 	import std.functional : binaryFun;
76 	import nxt.static_iota : iota;
77 	foreach (const i; iota!(0, indexes.length / 2)) /* TODO: use static foreach here instead? */
78 	{
79 		const j = indexes[2*i];
80 		const k = indexes[2*i + 1];
81 
82 		static assert(j >= 0, "First part of index pair " ~ i.stringof ~ " is negative");
83 		static assert(k >= 0, "Second part of index pair " ~ i.stringof ~ " is negative");
84 
85 		if (!binaryFun!less(r[j], r[k]))
86 			r.swapAt(j, k);
87 	}
88 }
89 
90 /** Largest length supported by network sort `networkSortUpTo`. */
91 private enum networkSortMaxLength = 22;
92 
93 /** Sort at most the first `n` elements of `r` using comparison `less` using
94  * a networking sort.
95  *
96  * Note: Sorting networks are not unique, not even optimal solutions.
97  *
98  * See_Also: http://stackoverflow.com/questions/3903086/standard-sorting-networks-for-small-values-of-n
99  * See_Also: http://www.cs.brandeis.edu/~hugues/sorting_networks.html
100  */
101 auto networkSortUpTo(uint n, alias less = "a < b", Range)(Range r)
102 if (n >= 2 &&
103 	isRandomAccessRange!Range)
104 in(r.length >= n) {
105 	auto s = r[0 .. n];
106 
107 	static if (n < 2)
108 		return; // no sorting needed
109 
110 	static	  if (n == 2)
111 		s.conditionalSwap!(less, Range, 0, 1);
112 	else static if (n == 3)
113 		s.conditionalSwap!(less, Range,
114 						   0,1,
115 						   1,2,
116 						   0,1);
117 	else static if (n == 4)
118 		s.conditionalSwap!(less, Range,
119 						   0,1, 2,3, // 2 in parallel
120 						   0,2, 1,3, // 2 in parallel
121 						   1,2);
122 	else static if (n == 5)
123 		s.conditionalSwap!(less, Range,
124 						   0,1, 3,4,  // 2 in parallel
125 						   0,2,
126 						   0,3, 1,2,  // 2 in parallel
127 						   2,3, 1,4,  // 2 in parallel
128 						   1,2, 3,4); // 2 in parallel
129 	else static if (n == 6)
130 		s.conditionalSwap!(less, Range,
131 						   0,1, 2,3, 4,5, // 3 in parallel
132 						   0,2, 1,4, 3,5, // 3 in parallel
133 						   0,1, 2,3, 4,5, // 3 in parallel
134 						   1,2, 3,4,	  // 2 in parallel
135 						   2,3);
136 	else static if (n == 8)	 // Bitonic Sorter. 6-steps: https://en.wikipedia.org/wiki/Bitonic_sorter
137 		s.conditionalSwap!(less, Range,
138 						   0,1, 2,3, 4,5, 6,7,
139 						   0,3, 4,7, 1,2, 5,6,
140 						   0,1, 2,3, 4,5, 6,7,
141 						   0,7, 1,6, 2,5, 3,4,
142 						   0,2, 1,3, 4,6, 5,7,
143 						   0,1, 2,3, 4,5, 6,7);
144 	else static if (n == 9)	 // R. W. Floyd.
145 		s.conditionalSwap!(less, Range,
146 						   0,1, 3,4, 6,7,
147 						   1,2, 4,5, 7,8,
148 						   0,1, 3,4, 6,7,
149 						   0,3,
150 						   3,6,
151 						   0,3, 1,4,
152 						   4,7,
153 						   1,4, 2,5,
154 						   5,8,
155 						   1,3, 2,5,
156 						   2,6, 5,7,
157 						   4,6,
158 						   2,4,
159 						   2,3, 5,6);
160 	else static if (n == 10)	// A. Waksman.
161 		s.conditionalSwap!(less, Range,
162 						   0,5, 1,6, 2,7, 3,8, 4,9,
163 						   0,3, 1,4, 5,8, 6,9,
164 						   0,2, 3,6, 7,9,
165 						   0,1, 2,4, 5,7, 8,9,
166 						   1,2, 3,5, 4,6, 7,8,
167 						   1,3, 2,5, 4,7, 6,8,
168 						   2,3, 6,7,
169 						   3,4, 5,6,
170 						   4,5);
171 	else static if (n == 11)	// 12-input by Shapiro and Green, minus the connections to a twelfth input.
172 		s.conditionalSwap!(less, Range,
173 						   0,1, 2,3, 4,5, 6,7, 8,9,
174 						   0,2, 1,3, 4,6, 5,7, 8,10,
175 						   1,2, 5,6, 9,10,
176 						   1,5, 6,10,
177 						   2,6, 5,9,
178 						   0,4, 1,5, 6,10,
179 						   3,7, 4,8,
180 						   0,4,
181 						   1,4, 3,8, 7,10,
182 						   2,3, 8,9,
183 						   2,4, 3,5, 6,8, 7,9,
184 						   3,4, 5,6, 7,8);
185 	else static if (n == 12)	// Shapiro and Green.
186 		s.conditionalSwap!(less, Range,
187 						   0,1, 2,3, 4,5, 6,7, 8,9, 10,11,
188 						   0,2, 1,3, 4,6, 5,7, 9,11, 8,10,
189 						   1,2, 5,6, 9,10,
190 						   1,5, 6,10,
191 						   2,6, 5,9,
192 						   0,4, 1,5, 6,10, 7,11,
193 						   3,7, 4,8,
194 						   0,4, 7,11,
195 						   1,4, 3,8, 7,10,
196 						   2,3, 8,9,
197 						   2,4, 3,5, 6,8, 7,9,
198 						   3,4, 5,6, 7,8);
199 	else static if (n == 13)	// Generated by the END algorithm.
200 		s.conditionalSwap!(less, Range,
201 						   0,12, 1,7, 2,6, 3,4, 5,8, 9,11,
202 						   0,1, 2,3, 4,6, 5,9, 7,12, 8,11,
203 						   0,2, 1,4, 3,7, 6,12, 7,8, 10,11,
204 						   4,9, 6,10, 11,12,
205 						   1,7, 3,4, 5,6, 8,9, 10,11,
206 						   1,3, 2,6, 4,7, 8,10, 9,11,
207 						   0,5, 2,5, 6,8, 9,10,
208 						   1,2, 3,5, 4,6, 7,8,
209 						   2,3, 4,5, 6,7, 8,9,
210 						   3,4, 5,6);
211 	else static if (n == 14) // Green's construction for 16 inputs minus connections to the fifteenth and sixteenth inputs.
212 		s.conditionalSwap!(less, Range,
213 						   0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13,
214 						   0,2, 1,3, 4,6, 5,7, 8,10, 9,11,
215 						   0,4, 1,5, 2,6, 3,7, 8,12, 9,13,
216 						   0,8, 1,9, 2,10, 3,11, 4,12, 5,13,
217 						   3,12, 5,10, 6,9,
218 						   1,2, 4,8, 7,11,
219 						   1,4, 2,8, 7,13,
220 						   2,4, 3,8, 5,6, 7,12, 9,10, 11,13,
221 						   3,5, 6,8, 7,9, 10,12,
222 						   3,4, 5,6, 7,8, 9,10, 11,12,
223 						   6,7, 8,9);
224 	else static if (n == 15) // Green's construction for 16 inputs minus connections to the sixteenth input.
225 		s.conditionalSwap!(less, Range,
226 						   0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13,
227 						   0,2, 1,3, 4,6, 5,7, 8,10, 9,11, 12,14,
228 						   0,4, 1,5, 2,6, 3,7, 8,12, 9,13, 10,14,
229 						   0,8, 1,9, 2,10, 3,11, 4,12, 5,13, 6,14,
230 						   3,12, 5,10, 6,9, 13,14,
231 						   1,2, 4,8, 7,11,
232 						   1,4, 2,8, 7,13, 11,14,
233 						   2,4, 3,8, 5,6, 7,12, 9,10, 11,13,
234 						   3,5, 6,8, 7,9, 10,12,
235 						   3,4, 5,6, 7,8, 9,10, 11,12,
236 						   6,7, 8,9);
237 	else static if (n == 16) // Green's construction. TODO: Use 10-step Bitonic sorter instead?
238 		s.conditionalSwap!(less, Range,
239 						   0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13, 14,15,
240 						   0,2, 1,3, 4,6, 5,7, 8,10, 9,11, 12,14, 13,15,
241 						   0,4, 1,5, 2,6, 3,7, 8,12, 9,13, 10,14, 11,15,
242 						   0,8, 1,9, 2,10, 3,11, 4,12, 5,13, 6,14, 7,15,
243 						   3,12, 5,10, 6,9, 13,14,
244 						   1,2, 4,8, 7,11,
245 						   1,4, 2,8, 7,13, 11,14,
246 						   2,4, 3,8, 5,6, 7,12, 9,10, 11,13,
247 						   3,5, 6,8, 7,9, 10,12,
248 						   3,4, 5,6, 7,8, 9,10, 11,12,
249 						   6,7, 8,9);
250 	else static if (n == 18) // Baddar's PHD thesis, chapter 6. Fewest stages but 2 comparators more than 'batcher'
251 		s.conditionalSwap!(less, Range,
252 						   0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13, 14,15, 16,17,
253 						   0,2, 1,3, 4,6, 5,7, 8,10, 9,11, 12,17, 13,14, 15,16,
254 						   0,4, 1,5, 2,6, 3,7, 9,10, 8,12, 11,16, 13,15, 14,17,
255 						   0,8, 1,13, 2,4, 3,5, 6,17, 7,16, 9,15, 10,14, 11,12,
256 						   0,1, 2,11, 3,15, 4,10, 5,12, 6,13, 7,14, 8,9, 16,17,
257 						   1,8, 3,10, 4,15, 5,11, 6,9, 7,13, 14,16,
258 						   1,2, 3,6, 4,8, 5,9, 7,11, 10,13, 12,16, 14,15,
259 						   2,3, 5,8, 6,7, 9,10, 11,14, 12,13,
260 						   2,4, 3,5, 6,8, 7,9, 10,11, 12,14, 13,15,
261 						   3,4, 5,6, 7,8, 9,10, 11,12, 13,14,
262 						   4,5, 6,7, 8,9, 10,11, 12,13);
263 	else static if (n == 22) // Baddar's PHD thesis, chapter 7. Fewest stages but 2 comparators more than 'batcher'
264 		s.conditionalSwap!(less, Range,
265 						   0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13, 14,15, 16,17, 18,19, 20,21,
266 						   0,5, 1,3, 2,4, 6,8, 7,9, 10,12, 11,13, 14,16, 15,17, 18,20, 19,21,
267 						   6,10, 7,11, 8,12, 9,13, 14,18, 15,19, 16,20, 17,21,
268 						   0,2, 1,4, 3,5, 7,15, 8,16, 9,17, 11,19,
269 						   0,10, 1,18, 3,12, 5,20, 13,21,
270 						   0,7, 2,4, 3,15, 6,14, 9,18, 17,20,
271 						   0,6, 1,8, 2,11, 3,8, 4,16, 5,10, 12,19, 13,14, 20,21,
272 						   2,13, 4,7, 5,9, 10,15, 11,17, 12,18, 14,16, 16,20, 18,19,
273 						   1,3, 2,6, 4,5, 7,9, 8,13, 10,11, 12,14, 15,17, 19,20,
274 						   3,5, 4,6, 7,8, 9,13, 10,12, 11,14, 15,18, 16,17,
275 						   1,2, 5,10, 6,7, 8,9, 11,12, 13,15, 14,16, 18,19,
276 						   2,3, 5,7, 8,10, 9,11, 12,13, 14,15, 16,18, 17,19,
277 						   2,4, 3,6, 7,8, 9,10, 11,12, 13,14, 15,16, 17,18,
278 						   3,4, 5,6, 8,9, 10,11, 12,13, 14,15, 16,17,
279 						   4,5, 6,7);
280 	else
281 		static assert(0, "Unsupported length: " ~ n.stringof);
282 
283 	import std.range : assumeSorted;
284 	return s.assumeSorted!less;
285 }
286 
287 /** Sort range `x` of length `n` using a networking sort.
288  */
289 auto networkSortExactly(uint n, alias less = "a < b", Range)(Range r)
290 if (n >= 2 &&
291 	isRandomAccessRange!Range)
292 in(r.length == n) {
293 	x[].networkSortUpTo!(n);
294 }
295 
296 /** Sort static array `x` of length `n` using a networking sort.
297  */
298 auto networkSortExactly(alias less = "a < b", T, size_t n)(ref T[n] x)
299 if (n >= 2) {
300 	x[].networkSortUpTo!(n);
301 }
302 
303 ///
304 pure nothrow @safe @nogc unittest {
305 	import std.algorithm.comparison : equal;
306 	int[4] x = [2, 3, 0, 1];
307 	const int[4] y = [0, 1, 2, 3];
308 	x.networkSortExactly();
309 	assert(x[].equal(y[]));
310 }
311 
312 /** Hybrid sort `r` using `networkSortUpTo` if length of `r` is
313  * less-than-or-equal to `networkSortMaxLength` and `std.algorithm.sorting.sort`
314  * otherwise.
315  *
316  * TODO: Add template parameter <= networkSortMaxLength that limits the amount
317  * networks generated.
318  */
319 auto hybridSort(alias less = "a < b", Range)(Range r)
320 if (isRandomAccessRange!Range) {
321 	import std.algorithm.sorting : isSorted;
322 	static foreach (n; 2 .. networkSortMaxLength + 1) {
323 		static if (__traits(compiles, { r.networkSortUpTo!(n, less); })) {
324 			if (n == r.length) {
325 				auto s = r.networkSortUpTo!(n, less);
326 				assert(s.isSorted!less);
327 				return s;
328 			}
329 		}
330 	}
331 	import std.algorithm.sorting : sort;
332 	return sort!less(r);
333 }
334 
335 ///
336 pure @safe unittest {
337 	import std.algorithm.sorting : isSorted;
338 	import std.algorithm.iteration : permutations;
339 	import std.range : iota;
340 	import std.random : randomShuffle, Random;
341 
342 	Random random;
343 
344 	alias T = uint;
345 
346 	const maxFullPermutationTestLength = 8;
347 	const maxTriedShufflings = 10_000; // maximum number of shufflings to try
348 
349 	import std.meta : AliasSeq;
350 	foreach (less; AliasSeq!("a < b", "a > b")) {
351 		foreach (const n; 0 .. networkSortMaxLength + 1) {
352 			if (n > maxFullPermutationTestLength) // if number of elements is too large
353 			{
354 				foreach (x; 0 .. maxTriedShufflings) {
355 					import std.array : array;
356 					auto y = iota(0, n).array;
357 					y.randomShuffle(random);
358 					y.hybridSort!less();
359 					assert(y.isSorted!less);
360 				}
361 			}
362 			else {
363 				foreach (x; iota(0, n).permutations) {
364 					import std.array : array;
365 					auto y = x.array;
366 					y.hybridSort!less();
367 					assert(y.isSorted!less);
368 				}
369 			}
370 		}
371 	}
372 }