libStatGen Software  1
Sort.cpp
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends