src/util/qsort_r.c

   1 /* Copyright (c) 2007-2014 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22 #include "nlopt-util.h"
  23 #include <stdlib.h>
  24
  25 #if defined(_MSC_VER)
  26 #define inline __inline
  27 #endif
  28
  29 typedef int              cmp_t(void *, const void *, const void *);
  30
  31 static inline char      *med3(char *, char *, char *, cmp_t *, void *);
  32
  33 #define MIN(a, b)       ((a) < (b) ? a : b)
  34
  35 /*
  36  * Qsort routine from Bentley & McIlroy's "Engineering a Sort Function".
  37  */
  38
  39 static inline void
  40 swapfunc(char *a, char *b, size_t es)
  41 {
  42         char t;
  43
  44         do {
  45                 t = *a;
  46                 *a++ = *b;
  47                 *b++ = t;
  48         } while (--es > 0);
  49 }
  50
  51 #define vecswap(a, b, n)                                \
  52         if ((n) > 0) swapfunc(a, b, n)
  53
  54 #define CMP(t, x, y) (cmp((t), (x), (y)))
  55
  56 static inline char *
  57 med3(char *a, char *b, char *c, cmp_t *cmp, void *thunk)
  58 {
  59         return CMP(thunk, a, b) < 0 ?
  60                (CMP(thunk, b, c) < 0 ? b : (CMP(thunk, a, c) < 0 ? c : a ))
  61               :(CMP(thunk, b, c) > 0 ? b : (CMP(thunk, a, c) < 0 ? a : c ));
  62 }
  63
  64
  65 void qsort_r_fallback(void *a, size_t n, size_t es, void *thunk, cmp_t *cmp)
  66 {
  67         char *pa, *pb, *pc, *pd, *pl, *pm, *pn;
  68         size_t d1, d2;
  69         int cmp_result;
  70         int swap_cnt;
  71
  72 loop:
  73         swap_cnt = 0;
  74         if (n < 7) {
  75                 for (pm = (char *)a + es; pm < (char *)a + n * es; pm += es)
  76                         for (pl = pm;
  77                              pl > (char *)a && CMP(thunk, pl - es, pl) > 0;
  78                              pl -= es)
  79                                 swapfunc(pl, pl - es, es);
  80                 return;
  81         }
  82         pm = (char *)a + (n / 2) * es;
  83         if (n > 7) {
  84                 pl = a;
  85                 pn = (char *)a + (n - 1) * es;
  86                 if (n > 40) {
  87                         size_t d = (n / 8) * es;
  88
  89                         pl = med3(pl, pl + d, pl + 2 * d, cmp, thunk);
  90                         pm = med3(pm - d, pm, pm + d, cmp, thunk);
  91                         pn = med3(pn - 2 * d, pn - d, pn, cmp, thunk);
  92                 }
  93                 pm = med3(pl, pm, pn, cmp, thunk);
  94         }
  95         swapfunc(a, pm, es);
  96         pa = pb = (char *)a + es;
  97
  98         pc = pd = (char *)a + (n - 1) * es;
  99         for (;;) {
 100                 while (pb <= pc && (cmp_result = CMP(thunk, pb, a)) <= 0) {
 101                         if (cmp_result == 0) {
 102                                 swap_cnt = 1;
 103                                 swapfunc(pa, pb, es);
 104                                 pa += es;
 105                         }
 106                         pb += es;
 107                 }
 108                 while (pb <= pc && (cmp_result = CMP(thunk, pc, a)) >= 0) {
 109                         if (cmp_result == 0) {
 110                                 swap_cnt = 1;
 111                                 swapfunc(pc, pd, es);
 112                                 pd -= es;
 113                         }
 114                         pc -= es;
 115                 }
 116                 if (pb > pc)
 117                         break;
 118                 swapfunc(pb, pc, es);
 119                 swap_cnt = 1;
 120                 pb += es;
 121                 pc -= es;
 122         }
 123         if (swap_cnt == 0) {  /* Switch to insertion sort */
 124                 for (pm = (char *)a + es; pm < (char *)a + n * es; pm += es)
 125                         for (pl = pm;
 126                              pl > (char *)a && CMP(thunk, pl - es, pl) > 0;
 127                              pl -= es)
 128                                 swapfunc(pl, pl - es, es);
 129                 return;
 130         }
 131
 132         pn = (char *)a + n * es;
 133         d1 = MIN(pa - (char *)a, pb - pa);
 134         vecswap(a, pb - d1, d1);
 135         d1 = MIN(pd - pc, pn - pd - es);
 136         vecswap(pb, pn - d1, d1);
 137
 138         d1 = pb - pa;
 139         d2 = pd - pc;
 140         if (d1 <= d2) {
 141                 /* Recurse on left partition, then iterate on right partition */
 142                 if (d1 > es) {
 143                         qsort_r_fallback(a, d1 / es, es, thunk, cmp);
 144                 }
 145                 if (d2 > es) {
 146                         /* Iterate rather than recurse to save stack space */
 147                         /* qsort(pn - d2, d2 / es, es, cmp); */
 148                         a = pn - d2;
 149                         n = d2 / es;
 150                         goto loop;
 151                 }
 152         } else {
 153                 /* Recurse on right partition, then iterate on left partition */
 154                 if (d2 > es) {
 155                         qsort_r_fallback(pn - d2, d2 / es, es, thunk, cmp);
 156                 }
 157                 if (d1 > es) {
 158                         /* Iterate rather than recurse to save stack space */
 159                         /* qsort(a, d1 / es, es, cmp); */
 160                         n = d1 / es;
 161                         goto loop;
 162                 }
 163         }
 164 }
 165
 166 /* these are required for GNU api compatibility as nlopt uses the BSD arguments ordering */
 167 typedef struct {
 168   cmp_t* compar;
 169   void *thunk;
 170 } qsort_wrapper;
 171
 172 static int qsort_cmp_wrap(const void *a, const void *b, void *thunk)
 173 {
 174   qsort_wrapper *wrap = (qsort_wrapper *) thunk;
 175   return (*wrap->compar)(wrap->thunk, a, b);
 176 }
 177
 178 void nlopt_qsort_r(void *base_, size_t nmemb, size_t size, void *thunk, cmp_t* compar)
 179 {
 180 #if defined(HAVE_QSORT_R) && (defined(__APPLE__) || defined(__FreeBSD__))
 181     qsort_r(base_, nmemb, size, thunk, compar);
 182 #elif defined(HAVE_QSORT_R) && defined(__linux__)
 183     qsort_wrapper wrapper;
 184     wrapper.compar = compar;
 185     wrapper.thunk = thunk;
 186     qsort_r(base_, nmemb, size, qsort_cmp_wrap, &wrapper);
 187 /*#elif defined(_WIN32)
 188   qsort_s(base_, nmemb, size, compar, thunk);*/
 189 #else
 190   qsort_r_fallback(base_, nmemb, size, thunk, compar);
 191 #endif
 192 }