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
34void 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
195void 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
346void * 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