lib/libc/stdlib/radixsort.c

   1 /*      $NetBSD: radixsort.c,v 1.18 2009/08/21 20:49:50 dsl Exp $       */
   2
   3 /*-
   4  * Copyright (c) 1990, 1993
   5  *      The Regents of the University of California.  All rights reserved.
   6  *
   7  * This code is derived from software contributed to Berkeley by
   8  * Peter McIlroy and by Dan Bernstein at New York University,
   9  *
  10  * Redistribution and use in source and binary forms, with or without
  11  * modification, are permitted provided that the following conditions
  12  * are met:
  13  * 1. Redistributions of source code must retain the above copyright
  14  *    notice, this list of conditions and the following disclaimer.
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in the
  17  *    documentation and/or other materials provided with the distribution.
  18  * 3. Neither the name of the University nor the names of its contributors
  19  *    may be used to endorse or promote products derived from this software
  20  *    without specific prior written permission.
  21  *
  22  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  23  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  24  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  25  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  26  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  27  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  28  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  29  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  30  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  31  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  32  * SUCH DAMAGE.
  33  */
  34
  35 #include <sys/cdefs.h>
  36 #if defined(LIBC_SCCS) && !defined(lint)
  37 #if 0
  38 static char sccsid[] = "@(#)radixsort.c 8.2 (Berkeley) 4/28/95";
  39 #else
  40 __RCSID("$NetBSD: radixsort.c,v 1.18 2009/08/21 20:49:50 dsl Exp $");
  41 #endif
  42 #endif /* LIBC_SCCS and not lint */
  43
  44 /*
  45  * Radixsort routines.
  46  *
  47  * Program r_sort_a() is unstable but uses O(logN) extra memory for a stack.
  48  * Use radixsort(a, n, trace, endchar) for this case.
  49  *
  50  * For stable sorting (using N extra pointers) use sradixsort(), which calls
  51  * r_sort_b().
  52  *
  53  * For a description of this code, see D. McIlroy, P. McIlroy, K. Bostic,
  54  * "Engineering Radix Sort".
  55  */
  56
  57 #include "namespace.h"
  58 #include <sys/types.h>
  59
  60 #include <assert.h>
  61 #include <errno.h>
  62 #include <stdlib.h>
  63
  64 #ifdef __weak_alias
  65 __weak_alias(radixsort,_radixsort)
  66 __weak_alias(sradixsort,_sradixsort)
  67 #endif
  68
  69 typedef struct {
  70         const u_char **sa;
  71         int sn, si;
  72 } stack;
  73
  74 static inline void simplesort(const u_char **, int, int, const u_char *, u_int);
  75 static void r_sort_a(const u_char **, int, int, const u_char *, u_int);
  76 static void r_sort_b(const u_char **,
  77             const u_char **, int, int, const u_char *, u_int);
  78
  79 #define THRESHOLD       20              /* Divert to simplesort(). */
  80 #define SIZE            512             /* Default stack size. */
  81
  82 #define SETUP {                                                         \
  83         if (tab == NULL) {                                              \
  84                 tr = tr0;                                               \
  85                 for (c = 0; c < endch; c++)                             \
  86                         tr0[c] = c + 1;                                 \
  87                 tr0[c] = 0;                                             \
  88                 for (c++; c < 256; c++)                                 \
  89                         tr0[c] = c;                                     \
  90                 endch = 0;                                              \
  91         } else {                                                        \
  92                 endch = tab[endch];                                     \
  93                 tr = tab;                                               \
  94                 if (endch != 0 && endch != 255) {                       \
  95                         errno = EINVAL;                                 \
  96                         return (-1);                                    \
  97                 }                                                       \
  98         }                                                               \
  99 }
 100
 101 int
 102 radixsort(const u_char **a, int n, const u_char *tab, u_int endch)
 103 {
 104         const u_char *tr;
 105         u_int c;
 106         u_char tr0[256];
 107
 108         _DIAGASSERT(a != NULL);
 109
 110         SETUP;
 111         r_sort_a(a, n, 0, tr, endch);
 112         return (0);
 113 }
 114
 115 int
 116 sradixsort(const u_char **a, int n, const u_char *tab, u_int endch)
 117 {
 118         const u_char *tr, **ta;
 119         u_int c;
 120         u_char tr0[256];
 121
 122         _DIAGASSERT(a != NULL);
 123         if (a == NULL) {
 124                 errno = EFAULT;
 125                 return (-1);
 126         }
 127
 128         SETUP;
 129         if (n < THRESHOLD)
 130                 simplesort(a, n, 0, tr, endch);
 131         else {
 132                 if ((ta = malloc(n * sizeof(a))) == NULL)
 133                         return (-1);
 134                 r_sort_b(a, ta, n, 0, tr, endch);
 135                 free(ta);
 136         }
 137         return (0);
 138 }
 139
 140 #define empty(s)        (s >= sp)
 141 #define pop(a, n, i)    a = (--sp)->sa, n = sp->sn, i = sp->si
 142 #define push(a, n, i)   sp->sa = a, sp->sn = n, (sp++)->si = i
 143 #define swap(a, b, t)   t = a, a = b, b = t
 144
 145 /* Unstable, in-place sort. */
 146 static void
 147 r_sort_a(const u_char **a, int n, int i, const u_char *tr, u_int endch)
 148 {
 149         static u_int count[256], nc, bmin;
 150         u_int c;
 151         const u_char **ak, *r;
 152         stack s[SIZE], *sp, *sp0, *sp1, temp;
 153         u_int *cp, bigc;
 154         const u_char **an, *t, **aj, **top[256];
 155
 156         _DIAGASSERT(a != NULL);
 157         _DIAGASSERT(tr != NULL);
 158
 159         /* Set up stack. */
 160         sp = s;
 161         push(a, n, i);
 162         while (!empty(s)) {
 163                 pop(a, n, i);
 164                 if (n < THRESHOLD) {
 165                         simplesort(a, n, i, tr, endch);
 166                         continue;
 167                 }
 168                 an = a + n;
 169
 170                 /* Make character histogram. */
 171                 if (nc == 0) {
 172                         bmin = 255;     /* First occupied bin, excluding eos. */
 173                         for (ak = a; ak < an;) {
 174                                 c = tr[(*ak++)[i]];
 175                                 if (++count[c] == 1 && c != endch) {
 176                                         if (c < bmin)
 177                                                 bmin = c;
 178                                         nc++;
 179                                 }
 180                         }
 181                         if (sp + nc > s + SIZE) {       /* Get more stack. */
 182                                 r_sort_a(a, n, i, tr, endch);
 183                                 continue;
 184                         }
 185                 }
 186
 187                 /*
 188                  * Set top[]; push incompletely sorted bins onto stack.
 189                  * top[] = pointers to last out-of-place element in bins.
 190                  * count[] = counts of elements in bins.
 191                  * Before permuting: top[c-1] + count[c] = top[c];
 192                  * during deal: top[c] counts down to top[c-1].
 193                  */
 194                 sp0 = sp1 = sp;         /* Stack position of biggest bin. */
 195                 bigc = 2;               /* Size of biggest bin. */
 196                 if (endch == 0)         /* Special case: set top[eos]. */
 197                         top[0] = ak = a + count[0];
 198                 else {
 199                         ak = a;
 200                         top[255] = an;
 201                 }
 202                 for (cp = count + bmin; nc > 0; cp++) {
 203                         while (*cp == 0)        /* Find next non-empty pile. */
 204                                 cp++;
 205                         if (*cp > 1) {
 206                                 if (*cp > bigc) {
 207                                         bigc = *cp;
 208                                         sp1 = sp;
 209                                 }
 210                                 push(ak, *cp, i+1);
 211                         }
 212                         top[cp-count] = ak += *cp;
 213                         nc--;
 214                 }
 215                 swap(*sp0, *sp1, temp); /* Play it safe -- biggest bin last. */
 216
 217                 /*
 218                  * Permute misplacements home.  Already home: everything
 219                  * before aj, and in bin[c], items from top[c] on.
 220                  * Inner loop:
 221                  *      r = next element to put in place;
 222                  *      ak = top[r[i]] = location to put the next element.
 223                  *      aj = bottom of 1st disordered bin.
 224                  * Outer loop:
 225                  *      Once the 1st disordered bin is done, ie. aj >= ak,
 226                  *      aj<-aj + count[c] connects the bins in a linked list;
 227                  *      reset count[c].
 228                  */
 229                 for (aj = a; aj < an;  *aj = r, aj += count[c], count[c] = 0)
 230                         for (r = *aj;  aj < (ak = --top[c = tr[r[i]]]);)
 231                                 swap(*ak, r, t);
 232         }
 233 }
 234
 235 /* Stable sort, requiring additional memory. */
 236 static void
 237 r_sort_b(const u_char **a, const u_char **ta, int n, int i, const u_char *tr,
 238     u_int endch)
 239 {
 240         static u_int count[256], nc, bmin;
 241         u_int c;
 242         const u_char **ak, **ai;
 243         stack s[512], *sp, *sp0, *sp1, temp;
 244         const u_char **top[256];
 245         u_int *cp, bigc;
 246
 247         _DIAGASSERT(a != NULL);
 248         _DIAGASSERT(ta != NULL);
 249         _DIAGASSERT(tr != NULL);
 250
 251         sp = s;
 252         push(a, n, i);
 253         while (!empty(s)) {
 254                 pop(a, n, i);
 255                 if (n < THRESHOLD) {
 256                         simplesort(a, n, i, tr, endch);
 257                         continue;
 258                 }
 259
 260                 if (nc == 0) {
 261                         bmin = 255;
 262                         for (ak = a + n; --ak >= a;) {
 263                                 c = tr[(*ak)[i]];
 264                                 if (++count[c] == 1 && c != endch) {
 265                                         if (c < bmin)
 266                                                 bmin = c;
 267                                         nc++;
 268                                 }
 269                         }
 270                         if (sp + nc > s + SIZE) {
 271                                 r_sort_b(a, ta, n, i, tr, endch);
 272                                 continue;
 273                         }
 274                 }
 275
 276                 sp0 = sp1 = sp;
 277                 bigc = 2;
 278                 if (endch == 0) {
 279                         top[0] = ak = a + count[0];
 280                         count[0] = 0;
 281                 } else {
 282                         ak = a;
 283                         top[255] = a + n;
 284                         count[255] = 0;
 285                 }
 286                 for (cp = count + bmin; nc > 0; cp++) {
 287                         while (*cp == 0)
 288                                 cp++;
 289                         if ((c = *cp) > 1) {
 290                                 if (c > bigc) {
 291                                         bigc = c;
 292                                         sp1 = sp;
 293                                 }
 294                                 push(ak, c, i+1);
 295                         }
 296                         top[cp-count] = ak += c;
 297                         *cp = 0;                        /* Reset count[]. */
 298                         nc--;
 299                 }
 300                 swap(*sp0, *sp1, temp);
 301
 302                 for (ak = ta + n, ai = a+n; ak > ta;)   /* Copy to temp. */
 303                         *--ak = *--ai;
 304                 for (ak = ta+n; --ak >= ta;)            /* Deal to piles. */
 305                         *--top[tr[(*ak)[i]]] = *ak;
 306         }
 307 }
 308
 309 /* insertion sort */
 310 static inline void
 311 simplesort(const u_char **a, int n, int b, const u_char *tr, u_int endch)
 312 {
 313         u_char ch;
 314         const u_char  **ak, **ai, *s, *t;
 315
 316         _DIAGASSERT(a != NULL);
 317         _DIAGASSERT(tr != NULL);
 318
 319         for (ak = a+1; --n >= 1; ak++)
 320                 for (ai = ak; ai > a; ai--) {
 321                         for (s = ai[0] + b, t = ai[-1] + b;
 322                             (ch = tr[*s]) != endch; s++, t++)
 323                                 if (ch != tr[*t])
 324                                         break;
 325                         if (ch >= tr[*t])
 326                                 break;
 327                         swap(ai[0], ai[-1], s);
 328                 }
 329 }