HermesCommon  3.0
qsort.cpp
Go to the documentation of this file.
1 // This file is part of HermesCommon.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
19 #include "qsort.h"
20 #include <limits.h>
21 #include <algorithm>
22 #include <stdlib.h>
23 
24 #define SWAP(a, b) { int c = *(a); *(a) = *(b); *(b) = c; }
25 
26 #define MAX_THRESH 4
27 
28 template<typename intType>
30 {
31 public:
32  intType *lo;
33  intType *hi;
34 };
35 
36 #define STACK_SIZE (CHAR_BIT * sizeof(size_t))
37 #define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
38 #define POP(low, high) ((void) (--top, (low = top->lo), (high = top->hi)))
39 #define STACK_NOT_EMPTY (stack < top)
40 
41 #define min(x, y) ((x) < (y) ? (x) : (y))
42 
43 template<typename intType>
44 void qsort_int(intType* pbase, size_t total_elems)
45 {
46  register intType *base_ptr = pbase;
47 
48  if (total_elems == 0) return;
49 
50  if (total_elems > MAX_THRESH)
51  {
52  intType *lo = base_ptr;
53  intType *hi = lo + total_elems - 1;
54  stack_node<intType> stack[STACK_SIZE];
55  stack_node<intType> *top = stack;
56 
57  PUSH(nullptr, nullptr);
58 
59  while (STACK_NOT_EMPTY)
60  {
61  intType *left_ptr;
62  intType *right_ptr;
63 
64  // Select median value from among LO, MID, and HI. Rearrange
65  // LO and HI so the three values are sorted. This lowers the
66  // probability of picking a pathological pivot value and
67  // skips a comparison for both the LEFT_PTR and RIGHT_PTR in
68  // the while loops.
69 
70  intType *mid = lo + ((hi - lo) >> 1);
71 
72  if (*mid < *lo)
73  SWAP(mid, lo);
74  if (*hi < *mid)
75  SWAP(mid, hi)
76  else
77  goto jump_over;
78  if (*mid < *lo)
79  SWAP(mid, lo);
80  jump_over:
81 
82  left_ptr = lo + 1;
83  right_ptr = hi - 1;
84 
85  // Here's the famous ``collapse the walls'' section of quicksort.
86  // Gotta like those tight inner loops! They are the main reason
87  // that this algorithm runs much faster than others.
88  do
89  {
90  while (*left_ptr < *mid)
91  left_ptr++;
92 
93  while (*mid < *right_ptr)
94  right_ptr--;
95 
96  if (left_ptr < right_ptr)
97  {
98  SWAP(left_ptr, right_ptr);
99  if (mid == left_ptr)
100  mid = right_ptr;
101  else if (mid == right_ptr)
102  mid = left_ptr;
103  left_ptr++;
104  right_ptr--;
105  }
106  else if (left_ptr == right_ptr)
107  {
108  left_ptr++;
109  right_ptr--;
110  break;
111  }
112  } while (left_ptr <= right_ptr);
113 
114  // Set up pointers for next iteration. First determine whether
115  // left and right partitions are below the threshold size. If so,
116  // ignore one or both. Otherwise, push the larger partition's
117  // bounds on the stack and continue sorting the smaller one.
118 
119  if ((size_t)(right_ptr - lo) <= MAX_THRESH)
120  {
121  if ((size_t)(hi - left_ptr) <= MAX_THRESH)
122  // Ignore both small partitions.
123  POP(lo, hi);
124  else
125  // Ignore small left partition.
126  lo = left_ptr;
127  }
128  else if ((size_t)(hi - left_ptr) <= MAX_THRESH)
129  // Ignore small right partition.
130  hi = right_ptr;
131  else if ((right_ptr - lo) > (hi - left_ptr))
132  {
133  // Push larger left partition indices.
134  PUSH(lo, right_ptr);
135  lo = left_ptr;
136  }
137  else
138  {
139  // Push larger right partition indices.
140  PUSH(left_ptr, hi);
141  hi = right_ptr;
142  }
143  }
144  }
145 
146  // Once the BASE_PTR array is partially sorted by quicksort the rest
147  // is completely sorted using insertion sort, since this is efficient
148  // for partitions below MAX_THRESH size. BASE_PTR points to the beginning
149  // of the array to sort, and END_PTR points at the very last element in
150  // the array (*not* one beyond it!).
151  {
152  intType *const end_ptr = base_ptr + total_elems - 1;
153  intType *tmp_ptr = base_ptr;
154  intType *thresh = min(end_ptr, base_ptr + MAX_THRESH);
155  register intType *run_ptr;
156 
157  // Find smallest element in first threshold and place it at the
158  // array's beginning. This is the smallest array element,
159  // and the operation speeds up insertion sort's inner loop.
160 
161  for (run_ptr = tmp_ptr + 1; run_ptr <= thresh; run_ptr++)
162  if (*run_ptr < *tmp_ptr)
163  tmp_ptr = run_ptr;
164 
165  if (tmp_ptr != base_ptr)
166  SWAP(tmp_ptr, base_ptr);
167 
168  // Insertion sort, running from left-hand-side up to right-hand-side.
169  run_ptr = base_ptr + 1;
170  while (++run_ptr <= end_ptr)
171  {
172  tmp_ptr = run_ptr - 1;
173  while (*run_ptr < *tmp_ptr)
174  tmp_ptr--;
175 
176  tmp_ptr++;
177  if (tmp_ptr != run_ptr)
178  {
179  intType c = *run_ptr;
180  register intType *hi, *lo;
181 
182  for (hi = lo = run_ptr; --lo >= tmp_ptr; hi = lo)
183  *hi = *lo;
184  *hi = c;
185  }
186  }
187  }
188 }
189 
190 template void qsort_int(int* pbase, size_t total_elems);
The QuickSort routine from glibc-2.5 modified for sorting int arrays.