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     {
63        void conditionalSwap(less = "a < b", Range)(Range r)
64        {
65        }
66     }
67 
68     instead.
69  */
70 private void conditionalSwap(alias less = "a < b", Range, indexes...)(Range r)
71 if (isRandomAccessRange!Range &&
72     allSatisfy!(isIntegral, typeof(indexes)) &&
73     indexes.length &&
74     (indexes.length & 1) == 0) // even number of indexes
75 {
76     import std.algorithm.mutation : swapAt;
77     import std.functional : binaryFun;
78     import nxt.static_iota : iota;
79     foreach (const i; iota!(0, indexes.length / 2)) /* TODO: use static foreach here instead? */
80     {
81         const j = indexes[2*i];
82         const k = indexes[2*i + 1];
83 
84         static assert(j >= 0, "First part of index pair " ~ i.stringof ~ " is negative");
85         static assert(k >= 0, "Second part of index pair " ~ i.stringof ~ " is negative");
86 
87         if (!binaryFun!less(r[j], r[k]))
88             r.swapAt(j, k);
89     }
90 }
91 
92 /** Largest length supported by network sort `networkSortUpTo`. */
93 private enum networkSortMaxLength = 22;
94 
95 /** Sort at most the first `n` elements of `r` using comparison `less` using
96  * a networking sort.
97  *
98  * Note: Sorting networks are not unique, not even optimal solutions.
99  *
100  * See_Also: http://stackoverflow.com/questions/3903086/standard-sorting-networks-for-small-values-of-n
101  * See_Also: http://www.cs.brandeis.edu/~hugues/sorting_networks.html
102  */
103 auto networkSortUpTo(uint n, alias less = "a < b", Range)(Range r)
104 if (n >= 2 &&
105     isRandomAccessRange!Range)
106 in(r.length >= n)
107 {
108     auto s = r[0 .. n];
109 
110     static if (n < 2)
111         return; // no sorting needed
112 
113     static      if (n == 2)
114         s.conditionalSwap!(less, Range, 0, 1);
115     else static if (n == 3)
116         s.conditionalSwap!(less, Range,
117                            0,1,
118                            1,2,
119                            0,1);
120     else static if (n == 4)
121         s.conditionalSwap!(less, Range,
122                            0,1, 2,3, // 2 in parallel
123                            0,2, 1,3, // 2 in parallel
124                            1,2);
125     else static if (n == 5)
126         s.conditionalSwap!(less, Range,
127                            0,1, 3,4,  // 2 in parallel
128                            0,2,
129                            0,3, 1,2,  // 2 in parallel
130                            2,3, 1,4,  // 2 in parallel
131                            1,2, 3,4); // 2 in parallel
132     else static if (n == 6)
133         s.conditionalSwap!(less, Range,
134                            0,1, 2,3, 4,5, // 3 in parallel
135                            0,2, 1,4, 3,5, // 3 in parallel
136                            0,1, 2,3, 4,5, // 3 in parallel
137                            1,2, 3,4,      // 2 in parallel
138                            2,3);
139     else static if (n == 8)     // Bitonic Sorter. 6-steps: https://en.wikipedia.org/wiki/Bitonic_sorter
140         s.conditionalSwap!(less, Range,
141                            0,1, 2,3, 4,5, 6,7,
142                            0,3, 4,7, 1,2, 5,6,
143                            0,1, 2,3, 4,5, 6,7,
144                            0,7, 1,6, 2,5, 3,4,
145                            0,2, 1,3, 4,6, 5,7,
146                            0,1, 2,3, 4,5, 6,7);
147     else static if (n == 9)     // R. W. Floyd.
148         s.conditionalSwap!(less, Range,
149                            0,1, 3,4, 6,7,
150                            1,2, 4,5, 7,8,
151                            0,1, 3,4, 6,7,
152                            0,3,
153                            3,6,
154                            0,3, 1,4,
155                            4,7,
156                            1,4, 2,5,
157                            5,8,
158                            1,3, 2,5,
159                            2,6, 5,7,
160                            4,6,
161                            2,4,
162                            2,3, 5,6);
163     else static if (n == 10)    // A. Waksman.
164         s.conditionalSwap!(less, Range,
165                            0,5, 1,6, 2,7, 3,8, 4,9,
166                            0,3, 1,4, 5,8, 6,9,
167                            0,2, 3,6, 7,9,
168                            0,1, 2,4, 5,7, 8,9,
169                            1,2, 3,5, 4,6, 7,8,
170                            1,3, 2,5, 4,7, 6,8,
171                            2,3, 6,7,
172                            3,4, 5,6,
173                            4,5);
174     else static if (n == 11)    // 12-input by Shapiro and Green, minus the connections to a twelfth input.
175         s.conditionalSwap!(less, Range,
176                            0,1, 2,3, 4,5, 6,7, 8,9,
177                            0,2, 1,3, 4,6, 5,7, 8,10,
178                            1,2, 5,6, 9,10,
179                            1,5, 6,10,
180                            2,6, 5,9,
181                            0,4, 1,5, 6,10,
182                            3,7, 4,8,
183                            0,4,
184                            1,4, 3,8, 7,10,
185                            2,3, 8,9,
186                            2,4, 3,5, 6,8, 7,9,
187                            3,4, 5,6, 7,8);
188     else static if (n == 12)    // Shapiro and Green.
189         s.conditionalSwap!(less, Range,
190                            0,1, 2,3, 4,5, 6,7, 8,9, 10,11,
191                            0,2, 1,3, 4,6, 5,7, 9,11, 8,10,
192                            1,2, 5,6, 9,10,
193                            1,5, 6,10,
194                            2,6, 5,9,
195                            0,4, 1,5, 6,10, 7,11,
196                            3,7, 4,8,
197                            0,4, 7,11,
198                            1,4, 3,8, 7,10,
199                            2,3, 8,9,
200                            2,4, 3,5, 6,8, 7,9,
201                            3,4, 5,6, 7,8);
202     else static if (n == 13)    // Generated by the END algorithm.
203         s.conditionalSwap!(less, Range,
204                            0,12, 1,7, 2,6, 3,4, 5,8, 9,11,
205                            0,1, 2,3, 4,6, 5,9, 7,12, 8,11,
206                            0,2, 1,4, 3,7, 6,12, 7,8, 10,11,
207                            4,9, 6,10, 11,12,
208                            1,7, 3,4, 5,6, 8,9, 10,11,
209                            1,3, 2,6, 4,7, 8,10, 9,11,
210                            0,5, 2,5, 6,8, 9,10,
211                            1,2, 3,5, 4,6, 7,8,
212                            2,3, 4,5, 6,7, 8,9,
213                            3,4, 5,6);
214     else static if (n == 14) // Green's construction for 16 inputs minus connections to the fifteenth and sixteenth inputs.
215         s.conditionalSwap!(less, Range,
216                            0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13,
217                            0,2, 1,3, 4,6, 5,7, 8,10, 9,11,
218                            0,4, 1,5, 2,6, 3,7, 8,12, 9,13,
219                            0,8, 1,9, 2,10, 3,11, 4,12, 5,13,
220                            3,12, 5,10, 6,9,
221                            1,2, 4,8, 7,11,
222                            1,4, 2,8, 7,13,
223                            2,4, 3,8, 5,6, 7,12, 9,10, 11,13,
224                            3,5, 6,8, 7,9, 10,12,
225                            3,4, 5,6, 7,8, 9,10, 11,12,
226                            6,7, 8,9);
227     else static if (n == 15) // Green's construction for 16 inputs minus connections to the sixteenth input.
228         s.conditionalSwap!(less, Range,
229                            0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13,
230                            0,2, 1,3, 4,6, 5,7, 8,10, 9,11, 12,14,
231                            0,4, 1,5, 2,6, 3,7, 8,12, 9,13, 10,14,
232                            0,8, 1,9, 2,10, 3,11, 4,12, 5,13, 6,14,
233                            3,12, 5,10, 6,9, 13,14,
234                            1,2, 4,8, 7,11,
235                            1,4, 2,8, 7,13, 11,14,
236                            2,4, 3,8, 5,6, 7,12, 9,10, 11,13,
237                            3,5, 6,8, 7,9, 10,12,
238                            3,4, 5,6, 7,8, 9,10, 11,12,
239                            6,7, 8,9);
240     else static if (n == 16) // Green's construction. TODO: Use 10-step Bitonic sorter instead?
241         s.conditionalSwap!(less, Range,
242                            0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13, 14,15,
243                            0,2, 1,3, 4,6, 5,7, 8,10, 9,11, 12,14, 13,15,
244                            0,4, 1,5, 2,6, 3,7, 8,12, 9,13, 10,14, 11,15,
245                            0,8, 1,9, 2,10, 3,11, 4,12, 5,13, 6,14, 7,15,
246                            3,12, 5,10, 6,9, 13,14,
247                            1,2, 4,8, 7,11,
248                            1,4, 2,8, 7,13, 11,14,
249                            2,4, 3,8, 5,6, 7,12, 9,10, 11,13,
250                            3,5, 6,8, 7,9, 10,12,
251                            3,4, 5,6, 7,8, 9,10, 11,12,
252                            6,7, 8,9);
253     else static if (n == 18) // Baddar's PHD thesis, chapter 6. Fewest stages but 2 comparators more than 'batcher'
254         s.conditionalSwap!(less, Range,
255                            0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13, 14,15, 16,17,
256                            0,2, 1,3, 4,6, 5,7, 8,10, 9,11, 12,17, 13,14, 15,16,
257                            0,4, 1,5, 2,6, 3,7, 9,10, 8,12, 11,16, 13,15, 14,17,
258                            0,8, 1,13, 2,4, 3,5, 6,17, 7,16, 9,15, 10,14, 11,12,
259                            0,1, 2,11, 3,15, 4,10, 5,12, 6,13, 7,14, 8,9, 16,17,
260                            1,8, 3,10, 4,15, 5,11, 6,9, 7,13, 14,16,
261                            1,2, 3,6, 4,8, 5,9, 7,11, 10,13, 12,16, 14,15,
262                            2,3, 5,8, 6,7, 9,10, 11,14, 12,13,
263                            2,4, 3,5, 6,8, 7,9, 10,11, 12,14, 13,15,
264                            3,4, 5,6, 7,8, 9,10, 11,12, 13,14,
265                            4,5, 6,7, 8,9, 10,11, 12,13);
266     else static if (n == 22) // Baddar's PHD thesis, chapter 7. Fewest stages but 2 comparators more than 'batcher'
267         s.conditionalSwap!(less, Range,
268                            0,1, 2,3, 4,5, 6,7, 8,9, 10,11, 12,13, 14,15, 16,17, 18,19, 20,21,
269                            0,5, 1,3, 2,4, 6,8, 7,9, 10,12, 11,13, 14,16, 15,17, 18,20, 19,21,
270                            6,10, 7,11, 8,12, 9,13, 14,18, 15,19, 16,20, 17,21,
271                            0,2, 1,4, 3,5, 7,15, 8,16, 9,17, 11,19,
272                            0,10, 1,18, 3,12, 5,20, 13,21,
273                            0,7, 2,4, 3,15, 6,14, 9,18, 17,20,
274                            0,6, 1,8, 2,11, 3,8, 4,16, 5,10, 12,19, 13,14, 20,21,
275                            2,13, 4,7, 5,9, 10,15, 11,17, 12,18, 14,16, 16,20, 18,19,
276                            1,3, 2,6, 4,5, 7,9, 8,13, 10,11, 12,14, 15,17, 19,20,
277                            3,5, 4,6, 7,8, 9,13, 10,12, 11,14, 15,18, 16,17,
278                            1,2, 5,10, 6,7, 8,9, 11,12, 13,15, 14,16, 18,19,
279                            2,3, 5,7, 8,10, 9,11, 12,13, 14,15, 16,18, 17,19,
280                            2,4, 3,6, 7,8, 9,10, 11,12, 13,14, 15,16, 17,18,
281                            3,4, 5,6, 8,9, 10,11, 12,13, 14,15, 16,17,
282                            4,5, 6,7);
283     else
284         static assert(0, "Unsupported length: " ~ n.stringof);
285 
286     import std.range : assumeSorted;
287     return s.assumeSorted!less;
288 }
289 
290 /** Sort range `x` of length `n` using a networking sort.
291  */
292 auto networkSortExactly(uint n, alias less = "a < b", Range)(Range r)
293 if (n >= 2 &&
294     isRandomAccessRange!Range)
295 in(r.length == n)
296 {
297     x[].networkSortUpTo!(n);
298 }
299 
300 /** Sort static array `x` of length `n` using a networking sort.
301  */
302 auto networkSortExactly(alias less = "a < b", T, size_t n)(ref T[n] x)
303 if (n >= 2)
304 {
305     x[].networkSortUpTo!(n);
306 }
307 
308 ///
309 @safe pure nothrow @nogc unittest
310 {
311     import std.algorithm.comparison : equal;
312     int[4] x = [2, 3, 0, 1];
313     const int[4] y = [0, 1, 2, 3];
314     x.networkSortExactly();
315     assert(x[].equal(y[]));
316 }
317 
318 /** Hybrid sort `r` using `networkSortUpTo` if length of `r` is
319  * less-than-or-equal to `networkSortMaxLength` and `std.algorithm.sorting.sort`
320  * otherwise.
321  *
322  * TODO: Add template parameter <= networkSortMaxLength that limits the amount
323  * networks generated.
324  */
325 auto hybridSort(alias less = "a < b", Range)(Range r)
326 if (isRandomAccessRange!Range)
327 {
328     import std.algorithm.sorting : isSorted;
329     static foreach (n; 2 .. networkSortMaxLength + 1)
330     {
331         static if (__traits(compiles, { r.networkSortUpTo!(n, less); }))
332         {
333             if (n == r.length)
334             {
335                 auto s = r.networkSortUpTo!(n, less);
336                 assert(s.isSorted!less);
337                 return s;
338             }
339         }
340     }
341     import std.algorithm.sorting : sort;
342     return sort!less(r);
343 }
344 
345 ///
346 @safe pure unittest
347 {
348     import std.algorithm.sorting : isSorted;
349     import std.algorithm.iteration : permutations;
350     import std.range : iota;
351     import std.random : randomShuffle, Random;
352 
353     Random random;
354 
355     alias T = uint;
356 
357     const maxFullPermutationTestLength = 8;
358     const maxTriedShufflings = 10_000; // maximum number of shufflings to try
359 
360     import std.meta : AliasSeq;
361     foreach (less; AliasSeq!("a < b", "a > b"))
362     {
363         foreach (const n; 0 .. networkSortMaxLength + 1)
364         {
365             if (n > maxFullPermutationTestLength) // if number of elements is too large
366             {
367                 foreach (x; 0 .. maxTriedShufflings)
368                 {
369                     import std.array : array;
370                     auto y = iota(0, n).array;
371                     y.randomShuffle(random);
372                     y.hybridSort!less();
373                     assert(y.isSorted!less);
374                 }
375             }
376             else
377             {
378                 foreach (x; iota(0, n).permutations)
379                 {
380                     import std.array : array;
381                     auto y = x.array;
382                     y.hybridSort!less();
383                     assert(y.isSorted!less);
384                 }
385             }
386         }
387     }
388 }