xref: /original-bsd/lib/libc/stdlib/radixsort.c (revision dfdcf295)
1 /*-
2  * Copyright (c) 1990 The Regents of the University of California.
3  * All rights reserved.
4  *
5  * %sccs.include.redist.c%
6  */
7 
8 #if defined(LIBC_SCCS) && !defined(lint)
9 static char sccsid[] = "@(#)radixsort.c	5.1 (Berkeley) 10/01/90";
10 #endif /* LIBC_SCCS and not lint */
11 
12 #include <sys/types.h>
13 #include <sys/errno.h>
14 #include <limits.h>
15 #include <stdlib.h>
16 #include <stddef.h>
17 
18 #define	NCHARS	(UCHAR_MAX + 1)
19 
20 /*
21  * shellsort (diminishing increment sort) from Data Structures and
22  * Algorithms, Aho, Hopcraft and Ullman, 1983 Edition, page 290;
23  * see also Knuth Vol. 3, page 84.  The increments are selected from
24  * formula (8), page 95.  Roughly O(N^3/2).
25  *
26  * __rspartition is the cutoff point for a further partitioning instead
27  * of a shellsort.  If it changes check __rsshell_increments.  Both of
28  * these are exported, as the best values are data dependent.  Unrolling
29  * this loop has not proven worthwhile.
30  */
31 #define	NPARTITION	40
32 int __rspartition = NPARTITION;
33 int __rsshell_increments[] = { 4, 1, 0, 0, 0, 0, 0, 0 };
34 #define SHELLSORT { \
35 	register u_char ch, *s1, *s2; \
36 	register int incr, *incrp; \
37 	for (incrp = __rsshell_increments; incr = *incrp++;) \
38 		for (t1 = incr; t1 < nmemb; ++t1) \
39 			for (t2 = t1 - incr; t2 >= 0;) { \
40 				s1 = p[t2] + indx; \
41 				s2 = p[t2 + incr] + indx; \
42 				while ((ch = tr[*s1++]) == tr[*s2] && ch) \
43 					++s2; \
44 				if (ch > tr[*s2]) { \
45 					s1 = p[t2]; \
46 					p[t2] = p[t2 + incr]; \
47 					p[t2 + incr] = s1; \
48 					t2 -= incr; \
49 				} else \
50 					break; \
51 			} \
52 }
53 
54 /*
55  * stack points to context structures.  Each structure defines a
56  * scheduled partitioning.  Radixsort exits when the stack is empty.
57  *
58  * The stack size is data dependent, and guessing is probably not
59  * worthwhile.  The initial stack fits in 1K with four bytes left over
60  * for malloc.  The initial size is exported, as the best value is
61  * data, and possibly, system, dependent.
62  */
63 typedef struct _stack {
64 	u_char **bot;
65 	int indx, nmemb;
66 } CONTEXT;
67 
68 int __radix_stacksize = (1024 - 4) / sizeof(CONTEXT);
69 #define	STACKPUSH { \
70 	if (stackp == estack) { \
71 		t1 = stackp - stack; \
72 		stackp = stack; \
73 		if (!(stack = (CONTEXT *)realloc((char *)stack, \
74 		    (__radix_stacksize *= 2) * sizeof(CONTEXT)))) { \
75 			t1 = errno; \
76 			free((char *)l2); \
77 			if (stackp) \
78 				free((char *)stackp); \
79 			errno = t1; \
80 			return(-1); \
81 		} \
82 		stackp = stack + t1; \
83 		estack = stack + __radix_stacksize; \
84 	} \
85 	stackp->bot = p; \
86 	stackp->nmemb = nmemb; \
87 	stackp->indx = indx; \
88 	++stackp; \
89 }
90 
91 #define	STACKPOP { \
92 	if (stackp == stack) \
93 		break; \
94 	--stackp; \
95 	bot = stackp->bot; \
96 	nmemb = stackp->nmemb; \
97 	indx = stackp->indx; \
98 }
99 
100 /*
101  * A variant of MSD radix sorting; see Knuth Vol. 3, page 177, and 5.2.5,
102  * Ex. 10 and 12.  Also, "Three Partition Refinement Algorithms, Paige and
103  * Tarjan, SIAM J. Comput. Vol. 16, No. 6, December 1987.
104  *
105  * This uses a simple sort as soon as a bucket crosses a cutoff point, rather
106  * than sorting the entire list after partitioning is finished.
107  *
108  * This is pure MSD instead of LSD of some number of MSD, switching to the
109  * simple sort as soon as possible.  Takes linear time relative to the number
110  * of bytes in the strings.
111  */
112 radixsort(l1, nmemb, tab, endbyte)
113 	u_char **l1, *tab, endbyte;
114 	register int nmemb;
115 {
116 	register int i, indx, t1, t2;
117 	register u_char **l2, **p, **bot, *tr;
118 	CONTEXT *estack, *stack, *stackp;
119 	int c[NCHARS + 1];
120 	u_char ltab[NCHARS];
121 
122 	if (nmemb <= 1)
123 		return(0);
124 
125 	/*
126 	 * there are two arrays, one provided by the user (l1), and the
127 	 * temporary one (l2).  The data is sorted to the temporary stack,
128 	 * and then copied back.  The speedup of using index to determine
129 	 * which stack the data is on and simply swapping stacks back and
130 	 * forth, thus avoiding the copy every iteration, turns out to not
131 	 * be any faster than the current implementation.
132 	 */
133 	if (!(l2 = (u_char **)malloc(sizeof(u_char *) * nmemb)))
134 		return(-1);
135 
136 	/* initialize stack */
137 	stack = stackp = estack = NULL;
138 
139 	/*
140 	 * tr references a table of sort weights; multiple entries may
141 	 * map to the same weight; EOS char must have the lowest weight.
142 	 */
143 	if (tab)
144 		tr = tab;
145 	else {
146 		tr = ltab;
147 		for (t1 = 0, t2 = endbyte; t1 < t2; ++t1)
148 			tr[t1] = t1 + 1;
149 		tr[t2] = 0;
150 		for (t1 = endbyte + 1; t1 < NCHARS; ++t1)
151 			tr[t1] = t1;
152 	}
153 
154 	/* first sort is entire stack */
155 	bot = l1;
156 	indx = 0;
157 
158 	for (;;) {
159 		/* clear bucket count array */
160 		bzero((char *)c, sizeof(c));
161 
162 		/*
163 		 * compute number of items that sort to the same bucket
164 		 * for this index.
165 		 */
166 		for (p = bot, i = nmemb; i--;)
167 			++c[tr[(*p++)[indx]]];
168 
169 		/*
170 		 * sum the number of characters into c, dividing the temp
171 		 * stack into the right number of buckets for this bucket,
172 		 * this index.  C contains the cumulative total of keys
173 		 * before and included in this bucket, and will later be
174 		 * used as an index to the bucket.  c[NCHARS] contains
175 		 * the total number of elements, for determining how many
176 		 * elements the last bucket contains.
177 		 */
178 		for (i = 1; i <= NCHARS; ++i)
179 			c[i] += c[i - 1];
180 
181 		/*
182 		 * partition the elements into buckets; c decrements
183 		 * through the bucket, and ends up pointing to the
184 		 * first element of the bucket.
185 		 */
186 		for (i = nmemb; i--;) {
187 			--p;
188 			l2[--c[tr[(*p)[indx]]]] = *p;
189 		}
190 
191 		/* copy the partitioned elements back to user stack */
192 		bcopy(l2, bot, nmemb * sizeof(u_char *));
193 
194 		++indx;
195 		/*
196 		 * sort buckets as necessary; don't sort c[0], it's the
197 		 * EOS character bucket, and nothing can follow EOS.
198 		 */
199 		for (i = NCHARS - 1; i; i--) {
200 			if ((nmemb = c[i + 1] - (t1 = c[i])) < 2)
201 				continue;
202 			p = bot + t1;
203 			if (nmemb > __rspartition)
204 				STACKPUSH
205 			else
206 				SHELLSORT
207 		}
208 		/* break out when stack is empty */
209 		STACKPOP
210 	}
211 
212 	free((char *)l2);
213 	free((char *)stack);
214 #ifdef STATS
215 	(void)fprintf(stderr, "max stack %u.\n", __radix_stacksize);
216 #endif
217 	return(0);
218 }
219