Home | History | Annotate | Line # | Download | only in stdlib
radixsort.c revision 1.16.30.1 
      1 /*	$NetBSD: radixsort.c,v 1.16.30.1 2009/05/13 19:18:27 jym 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.16.30.1 2009/05/13 19:18:27 jym 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
     75 	    __P((const u_char **, int, int, const u_char *, u_int));
     76 static void r_sort_a __P((const u_char **, int, int, const u_char *, u_int));
     77 static void r_sort_b __P((const u_char **,
     78 	    const u_char **, int, int, const u_char *, u_int));
     79 
     80 #define	THRESHOLD	20		/* Divert to simplesort(). */
     81 #define	SIZE		512		/* Default stack size. */
     82 
     83 #define SETUP {								\
     84 	if (tab == NULL) {						\
     85 		tr = tr0;						\
     86 		for (c = 0; c < endch; c++)				\
     87 			tr0[c] = c + 1;					\
     88 		tr0[c] = 0;						\
     89 		for (c++; c < 256; c++)					\
     90 			tr0[c] = c;					\
     91 		endch = 0;						\
     92 	} else {							\
     93 		endch = tab[endch];					\
     94 		tr = tab;						\
     95 		if (endch != 0 && endch != 255) {			\
     96 			errno = EINVAL;					\
     97 			return (-1);					\
     98 		}							\
     99 	}								\
    100 }
    101 
    102 int
    103 radixsort(a, n, tab, endch)
    104 	const u_char **a, *tab;
    105 	int n;
    106 	u_int endch;
    107 {
    108 	const u_char *tr;
    109 	u_int c;
    110 	u_char tr0[256];
    111 
    112 	_DIAGASSERT(a != NULL);
    113 	_DIAGASSERT(tab != NULL);
    114 
    115 	SETUP;
    116 	r_sort_a(a, n, 0, tr, endch);
    117 	return (0);
    118 }
    119 
    120 int
    121 sradixsort(a, n, tab, endch)
    122 	const u_char **a, *tab;
    123 	int n;
    124 	u_int endch;
    125 {
    126 	const u_char *tr, **ta;
    127 	u_int c;
    128 	u_char tr0[256];
    129 
    130 	_DIAGASSERT(a != NULL);
    131 	_DIAGASSERT(tab != NULL);
    132 	if (a == NULL || tab == NULL) {
    133 		errno = EFAULT;
    134 		return (-1);
    135 	}
    136 
    137 	SETUP;
    138 	if (n < THRESHOLD)
    139 		simplesort(a, n, 0, tr, endch);
    140 	else {
    141 		if ((ta = malloc(n * sizeof(a))) == NULL)
    142 			return (-1);
    143 		r_sort_b(a, ta, n, 0, tr, endch);
    144 		free(ta);
    145 	}
    146 	return (0);
    147 }
    148 
    149 #define empty(s)	(s >= sp)
    150 #define pop(a, n, i)	a = (--sp)->sa, n = sp->sn, i = sp->si
    151 #define push(a, n, i)	sp->sa = a, sp->sn = n, (sp++)->si = i
    152 #define swap(a, b, t)	t = a, a = b, b = t
    153 
    154 /* Unstable, in-place sort. */
    155 static void
    156 r_sort_a(a, n, i, tr, endch)
    157 	const u_char **a;
    158 	int n, i;
    159 	const u_char *tr;
    160 	u_int endch;
    161 {
    162 	static u_int count[256], nc, bmin;
    163 	u_int c;
    164 	const u_char **ak, *r;
    165 	stack s[SIZE], *sp, *sp0, *sp1, temp;
    166 	u_int *cp, bigc;
    167 	const u_char **an, *t, **aj, **top[256];
    168 
    169 	_DIAGASSERT(a != NULL);
    170 	_DIAGASSERT(tr != NULL);
    171 
    172 	/* Set up stack. */
    173 	sp = s;
    174 	push(a, n, i);
    175 	while (!empty(s)) {
    176 		pop(a, n, i);
    177 		if (n < THRESHOLD) {
    178 			simplesort(a, n, i, tr, endch);
    179 			continue;
    180 		}
    181 		an = a + n;
    182 
    183 		/* Make character histogram. */
    184 		if (nc == 0) {
    185 			bmin = 255;	/* First occupied bin, excluding eos. */
    186 			for (ak = a; ak < an;) {
    187 				c = tr[(*ak++)[i]];
    188 				if (++count[c] == 1 && c != endch) {
    189 					if (c < bmin)
    190 						bmin = c;
    191 					nc++;
    192 				}
    193 			}
    194 			if (sp + nc > s + SIZE) {	/* Get more stack. */
    195 				r_sort_a(a, n, i, tr, endch);
    196 				continue;
    197 			}
    198 		}
    199 
    200 		/*
    201 		 * Set top[]; push incompletely sorted bins onto stack.
    202 		 * top[] = pointers to last out-of-place element in bins.
    203 		 * count[] = counts of elements in bins.
    204 		 * Before permuting: top[c-1] + count[c] = top[c];
    205 		 * during deal: top[c] counts down to top[c-1].
    206 		 */
    207 		sp0 = sp1 = sp;		/* Stack position of biggest bin. */
    208 		bigc = 2;		/* Size of biggest bin. */
    209 		if (endch == 0)		/* Special case: set top[eos]. */
    210 			top[0] = ak = a + count[0];
    211 		else {
    212 			ak = a;
    213 			top[255] = an;
    214 		}
    215 		for (cp = count + bmin; nc > 0; cp++) {
    216 			while (*cp == 0)	/* Find next non-empty pile. */
    217 				cp++;
    218 			if (*cp > 1) {
    219 				if (*cp > bigc) {
    220 					bigc = *cp;
    221 					sp1 = sp;
    222 				}
    223 				push(ak, *cp, i+1);
    224 			}
    225 			top[cp-count] = ak += *cp;
    226 			nc--;
    227 		}
    228 		swap(*sp0, *sp1, temp);	/* Play it safe -- biggest bin last. */
    229 
    230 		/*
    231 		 * Permute misplacements home.  Already home: everything
    232 		 * before aj, and in bin[c], items from top[c] on.
    233 		 * Inner loop:
    234 		 *	r = next element to put in place;
    235 		 *	ak = top[r[i]] = location to put the next element.
    236 		 *	aj = bottom of 1st disordered bin.
    237 		 * Outer loop:
    238 		 *	Once the 1st disordered bin is done, ie. aj >= ak,
    239 		 *	aj<-aj + count[c] connects the bins in a linked list;
    240 		 *	reset count[c].
    241 		 */
    242 		for (aj = a; aj < an;  *aj = r, aj += count[c], count[c] = 0)
    243 			for (r = *aj;  aj < (ak = --top[c = tr[r[i]]]);)
    244 				swap(*ak, r, t);
    245 	}
    246 }
    247 
    248 /* Stable sort, requiring additional memory. */
    249 static void
    250 r_sort_b(a, ta, n, i, tr, endch)
    251 	const u_char **a, **ta;
    252 	int n, i;
    253 	const u_char *tr;
    254 	u_int endch;
    255 {
    256 	static u_int count[256], nc, bmin;
    257 	u_int c;
    258 	const u_char **ak, **ai;
    259 	stack s[512], *sp, *sp0, *sp1, temp;
    260 	const u_char **top[256];
    261 	u_int *cp, bigc;
    262 
    263 	_DIAGASSERT(a != NULL);
    264 	_DIAGASSERT(ta != NULL);
    265 	_DIAGASSERT(tr != NULL);
    266 
    267 	sp = s;
    268 	push(a, n, i);
    269 	while (!empty(s)) {
    270 		pop(a, n, i);
    271 		if (n < THRESHOLD) {
    272 			simplesort(a, n, i, tr, endch);
    273 			continue;
    274 		}
    275 
    276 		if (nc == 0) {
    277 			bmin = 255;
    278 			for (ak = a + n; --ak >= a;) {
    279 				c = tr[(*ak)[i]];
    280 				if (++count[c] == 1 && c != endch) {
    281 					if (c < bmin)
    282 						bmin = c;
    283 					nc++;
    284 				}
    285 			}
    286 			if (sp + nc > s + SIZE) {
    287 				r_sort_b(a, ta, n, i, tr, endch);
    288 				continue;
    289 			}
    290 		}
    291 
    292 		sp0 = sp1 = sp;
    293 		bigc = 2;
    294 		if (endch == 0) {
    295 			top[0] = ak = a + count[0];
    296 			count[0] = 0;
    297 		} else {
    298 			ak = a;
    299 			top[255] = a + n;
    300 			count[255] = 0;
    301 		}
    302 		for (cp = count + bmin; nc > 0; cp++) {
    303 			while (*cp == 0)
    304 				cp++;
    305 			if ((c = *cp) > 1) {
    306 				if (c > bigc) {
    307 					bigc = c;
    308 					sp1 = sp;
    309 				}
    310 				push(ak, c, i+1);
    311 			}
    312 			top[cp-count] = ak += c;
    313 			*cp = 0;			/* Reset count[]. */
    314 			nc--;
    315 		}
    316 		swap(*sp0, *sp1, temp);
    317 
    318 		for (ak = ta + n, ai = a+n; ak > ta;)	/* Copy to temp. */
    319 			*--ak = *--ai;
    320 		for (ak = ta+n; --ak >= ta;)		/* Deal to piles. */
    321 			*--top[tr[(*ak)[i]]] = *ak;
    322 	}
    323 }
    324 
    325 static inline void
    326 simplesort(a, n, b, tr, endch)	/* insertion sort */
    327 	const u_char **a;
    328 	int n, b;
    329 	const u_char *tr;
    330 	u_int endch;
    331 {
    332 	u_char ch;
    333 	const u_char  **ak, **ai, *s, *t;
    334 
    335 	_DIAGASSERT(a != NULL);
    336 	_DIAGASSERT(tr != NULL);
    337 
    338 	for (ak = a+1; --n >= 1; ak++)
    339 		for (ai = ak; ai > a; ai--) {
    340 			for (s = ai[0] + b, t = ai[-1] + b;
    341 			    (ch = tr[*s]) != endch; s++, t++)
    342 				if (ch != tr[*t])
    343 					break;
    344 			if (ch >= tr[*t])
    345 				break;
    346 			swap(ai[0], ai[-1], s);
    347 		}
    348 }
    349