xref: /csrg-svn/lib/libc/stdlib/radixsort.c (revision 45345)
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.2 (Berkeley) 10/12/90";
10 #endif /* LIBC_SCCS and not lint */
11 
12 #include <sys/types.h>
13 #include <limits.h>
14 #include <stdlib.h>
15 #include <stddef.h>
16 
17 #define	NCHARS	(UCHAR_MAX + 1)
18 
19 /*
20  * Shellsort (diminishing increment sort) from Data Structures and
21  * Algorithms, Aho, Hopcraft and Ullman, 1983 Edition, page 290;
22  * see also Knuth Vol. 3, page 84.  The increments are selected from
23  * formula (8), page 95.  Roughly O(N^3/2).
24  *
25  * __rspartition is the cutoff point for a further partitioning instead
26  * of a shellsort.  If it changes check __rsshell_increments.  Both of
27  * these are exported, as the best values are data dependent.  Unrolling
28  * this loop has not proven worthwhile.
29  */
30 #define	NPARTITION	40
31 int __rspartition = NPARTITION;
32 int __rsshell_increments[] = { 4, 1, 0, 0, 0, 0, 0, 0 };
33 #define SHELLSORT { \
34 	register u_char ch, *s1, *s2; \
35 	register int incr, *incrp; \
36 	for (incrp = __rsshell_increments; incr = *incrp++;) \
37 		for (t1 = incr; t1 < nmemb; ++t1) \
38 			for (t2 = t1 - incr; t2 >= 0;) { \
39 				s1 = p[t2] + indx; \
40 				s2 = p[t2 + incr] + indx; \
41 				while ((ch = tr[*s1++]) == tr[*s2] && ch) \
42 					++s2; \
43 				if (ch > tr[*s2]) { \
44 					s1 = p[t2]; \
45 					p[t2] = p[t2 + incr]; \
46 					p[t2 + incr] = s1; \
47 					t2 -= incr; \
48 				} else \
49 					break; \
50 			} \
51 }
52 
53 /*
54  * Stackp points to context structures, where each structure schedules a
55  * partitioning.  Radixsort exits when the stack is empty.
56  *
57  * If the buckets are placed on the stack randomly, the worst case is when:
58  *
59  *	(nbuckets - 1) contain (npartitions + 1) elements, with the last
60  *	bucket containing (nelements - ((npartitions + 1) * (nbuckets - 1))
61  *	keys.
62  *
63  * In this case, stack growth is bounded by:
64  *
65  *	(nelements / (npartitions + 1)) - 1
66  *
67  * Therefore, we force the largest bucket to be pushed on the stack first.
68  * Then the worst case is when:
69  *
70  * 	(nbuckets - 2) buckets contain (npartitions + 1) elements, with
71  *	the remaining elements split equally between the first bucket
72  *	pushed and the last bucket pushed.
73  *
74  * In this case, stack growth is bounded when:
75  *
76  *	for (partition_cnt = 0; nelements > npartitions; ++partition_cnt)
77  *		nelements =
78  *		    (nelements - (npartitions + 1) * (nbuckets - 2)) / 2;
79  * The bound is:
80  *
81  *	limit = partition_cnt * (nbuckets - 1);
82  */
83 typedef struct _stack {
84 	u_char **bot;
85 	int indx, nmemb;
86 } CONTEXT;
87 
88 #define	STACKPUSH { \
89 	stackp->bot = p; \
90 	stackp->nmemb = nmemb; \
91 	stackp->indx = indx; \
92 	++stackp; \
93 }
94 
95 #define	STACKPOP { \
96 	if (stackp == stack) \
97 		break; \
98 	--stackp; \
99 	bot = stackp->bot; \
100 	nmemb = stackp->nmemb; \
101 	indx = stackp->indx; \
102 }
103 
104 /*
105  * A variant of MSD radix sorting; see Knuth Vol. 3, page 177, and 5.2.5,
106  * Ex. 10 and 12.  Also, "Three Partition Refinement Algorithms, Paige
107  * and Tarjan, SIAM J. Comput. Vol. 16, No. 6, December 1987.
108  *
109  * This uses a simple sort as soon as a bucket crosses a cutoff point,
110  * rather than sorting the entire list after partitioning is finished.
111  * This should be an advantage.
112  *
113  * This is pure MSD instead of LSD of some number of MSD, switching to
114  * the simple sort as soon as possible.  Takes linear time relative to
115  * the number of bytes in the strings.
116  */
117 radixsort(l1, nmemb, tab, endbyte)
118 	u_char **l1, *tab, endbyte;
119 	register int nmemb;
120 {
121 	register int i, indx, t1, t2;
122 	register u_char **l2, **p, **bot, *tr;
123 	CONTEXT *stack, *stackp;
124 	int c[NCHARS + 1], max;
125 	u_char ltab[NCHARS];
126 
127 	if (nmemb <= 1)
128 		return(0);
129 
130 	/*
131 	 * T1 is the constant part of the equation, the number of elements
132 	 * represented on the stack between the top and bottom entries.
133 	 * Don't round as the divide by 2 rounds down (correct for value
134 	 * being subtracted).  The nelem value has to be rounded up before
135 	 * each divide because we want an upper bound.
136 	 */
137 	t1 = ((__rspartition + 1) * (UCHAR_MAX - 2)) >> 1;
138 	for (i = 0, t2 = nmemb; t2 > __rspartition; i += UCHAR_MAX - 1)
139 		t2 = (++t2 >> 1) - t1;
140 	if (i) {
141 		if (!(stack = stackp =
142 		    (CONTEXT *)malloc(i * sizeof(CONTEXT))))
143 			return(-1);
144 	} else
145 		stack = stackp = NULL;
146 
147 	/*
148 	 * There are two arrays, one provided by the user (l1), and the
149 	 * temporary one (l2).  The data is sorted to the temporary stack,
150 	 * and then copied back.  The speedup of using index to determine
151 	 * which stack the data is on and simply swapping stacks back and
152 	 * forth, thus avoiding the copy every iteration, turns out to not
153 	 * be any faster than the current implementation.
154 	 */
155 	if (!(l2 = (u_char **)malloc(sizeof(u_char *) * nmemb)))
156 		return(-1);
157 
158 	/*
159 	 * Tr references a table of sort weights; multiple entries may
160 	 * map to the same weight; EOS char must have the lowest weight.
161 	 */
162 	if (tab)
163 		tr = tab;
164 	else {
165 		tr = ltab;
166 		for (t1 = 0, t2 = endbyte; t1 < t2; ++t1)
167 			tr[t1] = t1 + 1;
168 		tr[t2] = 0;
169 		for (t1 = endbyte + 1; t1 < NCHARS; ++t1)
170 			tr[t1] = t1;
171 	}
172 
173 	/* First sort is entire stack */
174 	bot = l1;
175 	indx = 0;
176 
177 	for (;;) {
178 		/* Clear bucket count array */
179 		bzero((char *)c, sizeof(c));
180 
181 		/*
182 		 * Compute number of items that sort to the same bucket
183 		 * for this index.
184 		 */
185 		for (p = bot, i = nmemb; i--;)
186 			++c[tr[(*p++)[indx]]];
187 
188 		/*
189 		 * Sum the number of characters into c, dividing the temp
190 		 * stack into the right number of buckets for this bucket,
191 		 * this index.  C contains the cumulative total of keys
192 		 * before and included in this bucket, and will later be
193 		 * used as an index to the bucket.  c[NCHARS] contains
194 		 * the total number of elements, for determining how many
195 		 * elements the last bucket contains.  At the same time
196 		 * find the largest bucket so it gets handled first.
197 		 */
198 		for (i = 1, t2 = -1; i <= NCHARS; ++i) {
199 			if ((t1 = c[i - 1]) > t2) {
200 				t2 = t1;
201 				max = i;
202 			}
203 			c[i] += t1;
204 		}
205 
206 		/*
207 		 * Partition the elements into buckets; c decrements
208 		 * through the bucket, and ends up pointing to the
209 		 * first element of the bucket.
210 		 */
211 		for (i = nmemb; i--;) {
212 			--p;
213 			l2[--c[tr[(*p)[indx]]]] = *p;
214 		}
215 
216 		/* Copy the partitioned elements back to user stack */
217 		bcopy(l2, bot, nmemb * sizeof(u_char *));
218 
219 		++indx;
220 		/*
221 		 * Sort buckets as necessary; don't sort c[0], it's the
222 		 * EOS character bucket, and nothing can follow EOS.
223 		 */
224 		for (i = max; i; --i) {
225 			if ((nmemb = c[i + 1] - (t1 = c[i])) < 2)
226 				continue;
227 			p = bot + t1;
228 			if (nmemb > __rspartition)
229 				STACKPUSH
230 			else
231 				SHELLSORT
232 		}
233 		for (i = max + 1; i < NCHARS; ++i) {
234 			if ((nmemb = c[i + 1] - (t1 = c[i])) < 2)
235 				continue;
236 			p = bot + t1;
237 			if (nmemb > __rspartition)
238 				STACKPUSH
239 			else
240 				SHELLSORT
241 		}
242 		/* Break out when stack is empty */
243 		STACKPOP
244 	}
245 
246 	free((char *)l2);
247 	free((char *)stack);
248 	return(0);
249 }
250