libStatGen Software
1
|
00001 /* 00002 * Copyright (C) 2010 Regents of the University of Michigan 00003 * 00004 * This program is free software: you can redistribute it and/or modify 00005 * it under the terms of the GNU General Public License as published by 00006 * the Free Software Foundation, either version 3 of the License, or 00007 * (at your option) any later version. 00008 * 00009 * This program is distributed in the hope that it will be useful, 00010 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 * GNU General Public License for more details. 00013 * 00014 * You should have received a copy of the GNU General Public License 00015 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00016 */ 00017 00018 #include "Sort.h" 00019 #include "Error.h" 00020 00021 #include <stddef.h> 00022 #include <string.h> 00023 00024 00025 #define Item(b) (base_char+(b)*width) 00026 #define IsBefore(x,y) ((cmp(Item(x),Item(y)))<0) 00027 #define Exchange(x,y) {\ 00028 memcpy(tmp,Item(x),width);\ 00029 memcpy(Item(x),Item(y),width);\ 00030 memcpy(Item(y),tmp,width);\ 00031 } 00032 #define TRUE 1 00033 00034 void QuickSort(void *base, size_t nelem, size_t width, 00035 int (*cmp)(const void *, const void *)) 00036 { 00037 struct __QuickSortStack 00038 { 00039 size_t left, right; 00040 }; 00041 00042 if (nelem <= 1) 00043 return; 00044 00045 // Create a pseudo-stack to avoid recursion 00046 00047 char * base_char = (char *) base; 00048 const size_t stackSize = 128; 00049 00050 __QuickSortStack * stack = new __QuickSortStack[stackSize]; 00051 char * tmp = new char [width]; 00052 00053 if ((stack == NULL) || (tmp == NULL)) 00054 error("Out of memory in QuickSort routine"); 00055 00056 size_t stackIdx = 0; 00057 00058 // Size of minimum partition to median of three 00059 const size_t Threshold = 7; 00060 00061 // current partitions 00062 00063 size_t lsize, rsize; 00064 size_t l, mid, r; 00065 size_t scanl, scanr, pivot; 00066 00067 l = 0; 00068 r = nelem - 1; 00069 00070 while (TRUE) 00071 { 00072 while (r > l) 00073 { 00074 if (r - l > Threshold) 00075 // QuickSort : median of three partitioning 00076 { 00077 mid = (r + l) / 2; 00078 00079 // sort l, mid, and r 00080 if (IsBefore(mid, l)) 00081 Exchange(mid, l); 00082 00083 if (IsBefore(r, l)) 00084 Exchange(r, l); 00085 00086 if (IsBefore(r, mid)) 00087 Exchange(r, mid); 00088 00089 // set up for partitioning... 00090 pivot = r - 1; 00091 00092 Exchange(mid, pivot); 00093 00094 scanl = l + 1; 00095 scanr = r - 2; 00096 } 00097 else 00098 { 00099 // set up random partition -- faster 00100 pivot = r; 00101 scanl = l; 00102 scanr = r - 1; 00103 } 00104 00105 while (TRUE) 00106 { 00107 // scan from left for element >= pivot 00108 while ((scanl < r) && IsBefore(scanl, pivot)) 00109 ++scanl; 00110 00111 while ((scanr > l) && IsBefore(pivot, scanr)) 00112 --scanr; 00113 00114 // if scans have met, we are done 00115 if (scanl >= scanr) 00116 break; 00117 00118 Exchange(scanl, scanr); 00119 00120 if (scanl < r) 00121 ++scanl; 00122 00123 if (scanr > l) 00124 --scanr; 00125 } 00126 00127 // Exchange final element 00128 Exchange(pivot, scanl); 00129 00130 // Place largest partition on stack 00131 lsize = scanl - l; 00132 rsize = r - scanl; 00133 00134 if (lsize > rsize) 00135 { 00136 // if size is one we are done 00137 ++ stackIdx; 00138 00139 if (stackIdx == stackSize) 00140 error("Out of Stack in QuickSort routine"); 00141 00142 stack[stackIdx].left = l; 00143 stack[stackIdx].right = scanl - 1; 00144 00145 if (rsize != 0) 00146 l = scanl + 1; 00147 else 00148 break; 00149 } 00150 else 00151 { 00152 // if size is one we are done 00153 ++ stackIdx; 00154 00155 if (stackIdx == stackSize) 00156 error("Out of Stack in QuickSort routine"); 00157 00158 stack[stackIdx].left = scanl + 1; 00159 stack[stackIdx].right = r; 00160 00161 if (lsize != 0) 00162 r = scanl - 1; 00163 else 00164 break; 00165 } 00166 } 00167 00168 // iterate with values from stack 00169 if (stackIdx) 00170 { 00171 l = stack[stackIdx].left; 00172 r = stack[stackIdx].right; 00173 00174 --stackIdx; 00175 } 00176 else 00177 break; 00178 } 00179 00180 delete [] stack; 00181 delete [] tmp; 00182 } 00183 00184 #define Item2(b) (base_char2+(b)*width) 00185 #define Exchange2(x,y) {\ 00186 memcpy(tmp,Item(x),width);\ 00187 memcpy(Item(x),Item(y),width);\ 00188 memcpy(Item(y),tmp,width);\ 00189 memcpy(tmp,Item2(x),width);\ 00190 memcpy(Item2(x),Item2(y),width);\ 00191 memcpy(Item2(y),tmp,width);\ 00192 } 00193 00194 00195 void QuickSort2(void *base, void *base2, size_t nelem, size_t width, 00196 int (*cmp)(const void *, const void *)) 00197 { 00198 struct __QuickSortStack 00199 { 00200 size_t left, right; 00201 }; 00202 00203 if (nelem <= 1) 00204 return; 00205 00206 // Create a pseudo-stack to avoid recursion 00207 00208 char * base_char = (char *) base; 00209 char * base_char2 = (char *) base2; 00210 const size_t stackSize = 128; 00211 00212 __QuickSortStack * stack = new __QuickSortStack[stackSize]; 00213 char * tmp = new char [width]; 00214 00215 if ((stack == NULL) || (tmp == NULL)) 00216 error("Out of memory in QuickSort routine"); 00217 00218 size_t stackIdx = 0; 00219 00220 // Size of minimum partition to median of three 00221 const size_t Threshold = 7; 00222 00223 // current partitions 00224 00225 size_t lsize, rsize; 00226 size_t l, mid, r; 00227 size_t scanl, scanr, pivot; 00228 00229 l = 0; 00230 r = nelem - 1; 00231 00232 while (TRUE) 00233 { 00234 while (r > l) 00235 { 00236 if (r - l > Threshold) 00237 // QuickSort : median of three partitioning 00238 { 00239 mid = (r + l) / 2; 00240 00241 // sort l, mid, and r 00242 if (IsBefore(mid, l)) 00243 Exchange2(mid, l); 00244 00245 if (IsBefore(r, l)) 00246 Exchange2(r, l); 00247 00248 if (IsBefore(r, mid)) 00249 Exchange2(r, mid); 00250 00251 // set up for partitioning... 00252 pivot = r - 1; 00253 00254 Exchange2(mid, pivot); 00255 00256 scanl = l + 1; 00257 scanr = r - 2; 00258 } 00259 else 00260 { 00261 // set up random partition -- faster 00262 pivot = r; 00263 scanl = l; 00264 scanr = r - 1; 00265 } 00266 00267 while (TRUE) 00268 { 00269 // scan from left for element >= pivot 00270 while ((scanl < r) && IsBefore(scanl, pivot)) 00271 ++scanl; 00272 00273 while ((scanr > l) && IsBefore(pivot, scanr)) 00274 --scanr; 00275 00276 // if scans have met, we are done 00277 if (scanl >= scanr) 00278 break; 00279 00280 Exchange2(scanl, scanr); 00281 00282 if (scanl < r) 00283 ++scanl; 00284 00285 if (scanr > l) 00286 --scanr; 00287 } 00288 00289 // Exchange final element 00290 Exchange2(pivot, scanl); 00291 00292 // Place largest partition on stack 00293 lsize = scanl - l; 00294 rsize = r - scanl; 00295 00296 if (lsize > rsize) 00297 { 00298 // if size is one we are done 00299 ++ stackIdx; 00300 00301 if (stackIdx == stackSize) 00302 error("Out of Stack in QuickSort routine"); 00303 00304 stack[stackIdx].left = l; 00305 stack[stackIdx].right = scanl - 1; 00306 00307 if (rsize != 0) 00308 l = scanl + 1; 00309 else 00310 break; 00311 } 00312 else 00313 { 00314 // if size is one we are done 00315 ++ stackIdx; 00316 00317 if (stackIdx == stackSize) 00318 error("Out of Stack in QuickSort routine"); 00319 00320 stack[stackIdx].left = scanl + 1; 00321 stack[stackIdx].right = r; 00322 00323 if (lsize != 0) 00324 r = scanl - 1; 00325 else 00326 break; 00327 } 00328 } 00329 00330 // iterate with values from stack 00331 if (stackIdx) 00332 { 00333 l = stack[stackIdx].left; 00334 r = stack[stackIdx].right; 00335 00336 --stackIdx; 00337 } 00338 else 00339 break; 00340 } 00341 00342 delete [] stack; 00343 delete [] tmp; 00344 } 00345 00346 void * BinarySearch(const void *key, const void *base, 00347 size_t nelem, size_t width, 00348 int (*cmp)(const void *, const void *)) 00349 { 00350 if (nelem == 0) 00351 return NULL; 00352 00353 char * base_char = (char *) base; 00354 00355 int left = 0; 00356 int right = nelem - 1; 00357 00358 while (right >= left) 00359 { 00360 int probe = (left + right) / 2; 00361 int test = cmp(key, Item(probe)); 00362 00363 if (test == 0) 00364 return (void *) Item(probe); 00365 00366 if (test < 0) 00367 right = probe - 1; 00368 else 00369 left = probe + 1; 00370 } 00371 00372 return NULL; 00373 } 00374