libStatGen Software  1
Sort.cpp
1 /*
2  * Copyright (C) 2010 Regents of the University of Michigan
3  *
4  * This program is free software: you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation, either version 3 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program. If not, see <http://www.gnu.org/licenses/>.
16  */
17 
18 #include "Sort.h"
19 #include "Error.h"
20 
21 #include <stddef.h>
22 #include <string.h>
23 
24 
25 #define Item(b) (base_char+(b)*width)
26 #define IsBefore(x,y) ((cmp(Item(x),Item(y)))<0)
27 #define Exchange(x,y) {\
28  memcpy(tmp,Item(x),width);\
29  memcpy(Item(x),Item(y),width);\
30  memcpy(Item(y),tmp,width);\
31  }
32 #define TRUE 1
33 
34 void QuickSort(void *base, size_t nelem, size_t width,
35  int (*cmp)(const void *, const void *))
36 {
37  struct __QuickSortStack
38  {
39  size_t left, right;
40  };
41 
42  if (nelem <= 1)
43  return;
44 
45  // Create a pseudo-stack to avoid recursion
46 
47  char * base_char = (char *) base;
48  const size_t stackSize = 128;
49 
50  __QuickSortStack * stack = new __QuickSortStack[stackSize];
51  char * tmp = new char [width];
52 
53  if ((stack == NULL) || (tmp == NULL))
54  error("Out of memory in QuickSort routine");
55 
56  size_t stackIdx = 0;
57 
58  // Size of minimum partition to median of three
59  const size_t Threshold = 7;
60 
61  // current partitions
62 
63  size_t lsize, rsize;
64  size_t l, mid, r;
65  size_t scanl, scanr, pivot;
66 
67  l = 0;
68  r = nelem - 1;
69 
70  while (TRUE)
71  {
72  while (r > l)
73  {
74  if (r - l > Threshold)
75  // QuickSort : median of three partitioning
76  {
77  mid = (r + l) / 2;
78 
79  // sort l, mid, and r
80  if (IsBefore(mid, l))
81  Exchange(mid, l);
82 
83  if (IsBefore(r, l))
84  Exchange(r, l);
85 
86  if (IsBefore(r, mid))
87  Exchange(r, mid);
88 
89  // set up for partitioning...
90  pivot = r - 1;
91 
92  Exchange(mid, pivot);
93 
94  scanl = l + 1;
95  scanr = r - 2;
96  }
97  else
98  {
99  // set up random partition -- faster
100  pivot = r;
101  scanl = l;
102  scanr = r - 1;
103  }
104 
105  while (TRUE)
106  {
107  // scan from left for element >= pivot
108  while ((scanl < r) && IsBefore(scanl, pivot))
109  ++scanl;
110 
111  while ((scanr > l) && IsBefore(pivot, scanr))
112  --scanr;
113 
114  // if scans have met, we are done
115  if (scanl >= scanr)
116  break;
117 
118  Exchange(scanl, scanr);
119 
120  if (scanl < r)
121  ++scanl;
122 
123  if (scanr > l)
124  --scanr;
125  }
126 
127  // Exchange final element
128  Exchange(pivot, scanl);
129 
130  // Place largest partition on stack
131  lsize = scanl - l;
132  rsize = r - scanl;
133 
134  if (lsize > rsize)
135  {
136  // if size is one we are done
137  ++ stackIdx;
138 
139  if (stackIdx == stackSize)
140  error("Out of Stack in QuickSort routine");
141 
142  stack[stackIdx].left = l;
143  stack[stackIdx].right = scanl - 1;
144 
145  if (rsize != 0)
146  l = scanl + 1;
147  else
148  break;
149  }
150  else
151  {
152  // if size is one we are done
153  ++ stackIdx;
154 
155  if (stackIdx == stackSize)
156  error("Out of Stack in QuickSort routine");
157 
158  stack[stackIdx].left = scanl + 1;
159  stack[stackIdx].right = r;
160 
161  if (lsize != 0)
162  r = scanl - 1;
163  else
164  break;
165  }
166  }
167 
168  // iterate with values from stack
169  if (stackIdx)
170  {
171  l = stack[stackIdx].left;
172  r = stack[stackIdx].right;
173 
174  --stackIdx;
175  }
176  else
177  break;
178  }
179 
180  delete [] stack;
181  delete [] tmp;
182 }
183 
184 #define Item2(b) (base_char2+(b)*width)
185 #define Exchange2(x,y) {\
186  memcpy(tmp,Item(x),width);\
187  memcpy(Item(x),Item(y),width);\
188  memcpy(Item(y),tmp,width);\
189  memcpy(tmp,Item2(x),width);\
190  memcpy(Item2(x),Item2(y),width);\
191  memcpy(Item2(y),tmp,width);\
192  }
193 
194 
195 void QuickSort2(void *base, void *base2, size_t nelem, size_t width,
196  int (*cmp)(const void *, const void *))
197 {
198  struct __QuickSortStack
199  {
200  size_t left, right;
201  };
202 
203  if (nelem <= 1)
204  return;
205 
206  // Create a pseudo-stack to avoid recursion
207 
208  char * base_char = (char *) base;
209  char * base_char2 = (char *) base2;
210  const size_t stackSize = 128;
211 
212  __QuickSortStack * stack = new __QuickSortStack[stackSize];
213  char * tmp = new char [width];
214 
215  if ((stack == NULL) || (tmp == NULL))
216  error("Out of memory in QuickSort routine");
217 
218  size_t stackIdx = 0;
219 
220  // Size of minimum partition to median of three
221  const size_t Threshold = 7;
222 
223  // current partitions
224 
225  size_t lsize, rsize;
226  size_t l, mid, r;
227  size_t scanl, scanr, pivot;
228 
229  l = 0;
230  r = nelem - 1;
231 
232  while (TRUE)
233  {
234  while (r > l)
235  {
236  if (r - l > Threshold)
237  // QuickSort : median of three partitioning
238  {
239  mid = (r + l) / 2;
240 
241  // sort l, mid, and r
242  if (IsBefore(mid, l))
243  Exchange2(mid, l);
244 
245  if (IsBefore(r, l))
246  Exchange2(r, l);
247 
248  if (IsBefore(r, mid))
249  Exchange2(r, mid);
250 
251  // set up for partitioning...
252  pivot = r - 1;
253 
254  Exchange2(mid, pivot);
255 
256  scanl = l + 1;
257  scanr = r - 2;
258  }
259  else
260  {
261  // set up random partition -- faster
262  pivot = r;
263  scanl = l;
264  scanr = r - 1;
265  }
266 
267  while (TRUE)
268  {
269  // scan from left for element >= pivot
270  while ((scanl < r) && IsBefore(scanl, pivot))
271  ++scanl;
272 
273  while ((scanr > l) && IsBefore(pivot, scanr))
274  --scanr;
275 
276  // if scans have met, we are done
277  if (scanl >= scanr)
278  break;
279 
280  Exchange2(scanl, scanr);
281 
282  if (scanl < r)
283  ++scanl;
284 
285  if (scanr > l)
286  --scanr;
287  }
288 
289  // Exchange final element
290  Exchange2(pivot, scanl);
291 
292  // Place largest partition on stack
293  lsize = scanl - l;
294  rsize = r - scanl;
295 
296  if (lsize > rsize)
297  {
298  // if size is one we are done
299  ++ stackIdx;
300 
301  if (stackIdx == stackSize)
302  error("Out of Stack in QuickSort routine");
303 
304  stack[stackIdx].left = l;
305  stack[stackIdx].right = scanl - 1;
306 
307  if (rsize != 0)
308  l = scanl + 1;
309  else
310  break;
311  }
312  else
313  {
314  // if size is one we are done
315  ++ stackIdx;
316 
317  if (stackIdx == stackSize)
318  error("Out of Stack in QuickSort routine");
319 
320  stack[stackIdx].left = scanl + 1;
321  stack[stackIdx].right = r;
322 
323  if (lsize != 0)
324  r = scanl - 1;
325  else
326  break;
327  }
328  }
329 
330  // iterate with values from stack
331  if (stackIdx)
332  {
333  l = stack[stackIdx].left;
334  r = stack[stackIdx].right;
335 
336  --stackIdx;
337  }
338  else
339  break;
340  }
341 
342  delete [] stack;
343  delete [] tmp;
344 }
345 
346 void * BinarySearch(const void *key, const void *base,
347  size_t nelem, size_t width,
348  int (*cmp)(const void *, const void *))
349 {
350  if (nelem == 0)
351  return NULL;
352 
353  char * base_char = (char *) base;
354 
355  int left = 0;
356  int right = nelem - 1;
357 
358  while (right >= left)
359  {
360  int probe = (left + right) / 2;
361  int test = cmp(key, Item(probe));
362 
363  if (test == 0)
364  return (void *) Item(probe);
365 
366  if (test < 0)
367  right = probe - 1;
368  else
369  left = probe + 1;
370  }
371 
372  return NULL;
373 }
374