[unix-history] / lib / libc / stdlib / radixsort.c

/*-
 * Copyright (c) 1990 The Regents of the University of California.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *	This product includes software developed by the University of
 *	California, Berkeley and its contributors.
 * 4. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#if defined(LIBC_SCCS) && !defined(lint)
static char sccsid[] = "@(#)radixsort.c	5.7 (Berkeley) 2/23/91";
#endif /* LIBC_SCCS and not lint */

#include <sys/types.h>
#include <limits.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>

/*
 * __rspartition is the cutoff point for a further partitioning instead
 * of a shellsort.  If it changes check __rsshell_increments.  Both of
 * these are exported, as the best values are data dependent.
 */
#define	NPARTITION	40
int __rspartition = NPARTITION;
int __rsshell_increments[] = { 4, 1, 0, 0, 0, 0, 0, 0 };

/*
 * Stackp points to context structures, where each structure schedules a
 * partitioning.  Radixsort exits when the stack is empty.
 *
 * If the buckets are placed on the stack randomly, the worst case is when
 * all the buckets but one contain (npartitions + 1) elements and the bucket
 * pushed on the stack last contains the rest of the elements.  In this case,
 * stack growth is bounded by:
 *
 *	limit = (nelements / (npartitions + 1)) - 1;
 *
 * This is a very large number, 52,377,648 for the maximum 32-bit signed int.
 *
 * By forcing the largest bucket to be pushed on the stack first, the worst
 * case is when all but two buckets each contain (npartitions + 1) elements,
 * with the remaining elements split equally between the first and last
 * buckets pushed on the stack.  In this case, stack growth is bounded when:
 *
 *	for (partition_cnt = 0; nelements > npartitions; ++partition_cnt)
 *		nelements =
 *		    (nelements - (npartitions + 1) * (nbuckets - 2)) / 2;
 * The bound is:
 *
 *	limit = partition_cnt * (nbuckets - 1);
 *
 * This is a much smaller number, 4590 for the maximum 32-bit signed int.
 */
#define	NBUCKETS	(UCHAR_MAX + 1)

typedef struct _stack {
	const u_char **bot;
	int indx, nmemb;
} CONTEXT;

#define	STACKPUSH { \
	stackp->bot = p; \
	stackp->nmemb = nmemb; \
	stackp->indx = indx; \
	++stackp; \
}
#define	STACKPOP { \
	if (stackp == stack) \
		break; \
	--stackp; \
	bot = stackp->bot; \
	nmemb = stackp->nmemb; \
	indx = stackp->indx; \
}

/*
 * A variant of MSD radix sorting; see Knuth Vol. 3, page 177, and 5.2.5,
 * Ex. 10 and 12.  Also, "Three Partition Refinement Algorithms, Paige
 * and Tarjan, SIAM J. Comput. Vol. 16, No. 6, December 1987.
 *
 * This uses a simple sort as soon as a bucket crosses a cutoff point,
 * rather than sorting the entire list after partitioning is finished.
 * This should be an advantage.
 *
 * This is pure MSD instead of LSD of some number of MSD, switching to
 * the simple sort as soon as possible.  Takes linear time relative to
 * the number of bytes in the strings.
 */
int
#if __STDC__
radixsort(const u_char **l1, int nmemb, const u_char *tab, u_char endbyte)
#else
radixsort(l1, nmemb, tab, endbyte)
	const u_char **l1;
	register int nmemb;
	const u_char *tab;
	u_char endbyte;
#endif
{
	register int i, indx, t1, t2;
	register const u_char **l2;
	register const u_char **p;
	register const u_char **bot;
	register const u_char *tr;
	CONTEXT *stack, *stackp;
	int c[NBUCKETS + 1], max;
	u_char ltab[NBUCKETS];
	static void shellsort();

	if (nmemb <= 1)
		return(0);

	/*
	 * T1 is the constant part of the equation, the number of elements
	 * represented on the stack between the top and bottom entries.
	 * It doesn't get rounded as the divide by 2 rounds down (correct
	 * for a value being subtracted).  T2, the nelem value, has to be
	 * rounded up before each divide because we want an upper bound;
	 * this could overflow if nmemb is the maximum int.
	 */
	t1 = ((__rspartition + 1) * (NBUCKETS - 2)) >> 1;
	for (i = 0, t2 = nmemb; t2 > __rspartition; i += NBUCKETS - 1)
		t2 = ((t2 + 1) >> 1) - t1;
	if (i) {
		if (!(stack = stackp = (CONTEXT *)malloc(i * sizeof(CONTEXT))))
			return(-1);
	} else
		stack = stackp = NULL;

	/*
	 * There are two arrays, one provided by the user (l1), and the
	 * temporary one (l2).  The data is sorted to the temporary stack,
	 * and then copied back.  The speedup of using index to determine
	 * which stack the data is on and simply swapping stacks back and
	 * forth, thus avoiding the copy every iteration, turns out to not
	 * be any faster than the current implementation.
	 */
	if (!(l2 = (const u_char **)malloc(sizeof(u_char *) * nmemb)))
		return(-1);

	/*
	 * Tr references a table of sort weights; multiple entries may
	 * map to the same weight; EOS char must have the lowest weight.
	 */
	if (tab)
		tr = tab;
	else {
		for (t1 = 0, t2 = endbyte; t1 < t2; ++t1)
			ltab[t1] = t1 + 1;
		ltab[t2] = 0;
		for (t1 = endbyte + 1; t1 < NBUCKETS; ++t1)
			ltab[t1] = t1;
		tr = ltab;
	}

	/* First sort is entire stack */
	bot = l1;
	indx = 0;

	for (;;) {
		/* Clear bucket count array */
		bzero((char *)c, sizeof(c));

		/*
		 * Compute number of items that sort to the same bucket
		 * for this index.
		 */
		for (p = bot, i = nmemb; --i >= 0;)
			++c[tr[(*p++)[indx]]];

		/*
		 * Sum the number of characters into c, dividing the temp
		 * stack into the right number of buckets for this bucket,
		 * this index.  C contains the cumulative total of keys
		 * before and included in this bucket, and will later be
		 * used as an index to the bucket.  c[NBUCKETS] contains
		 * the total number of elements, for determining how many
		 * elements the last bucket contains.  At the same time
		 * find the largest bucket so it gets pushed first.
		 */
		for (i = max = t1 = 0, t2 = __rspartition; i <= NBUCKETS; ++i) {
			if (c[i] > t2) {
				t2 = c[i];
				max = i;
			}
			t1 = c[i] += t1;
		}

		/*
		 * Partition the elements into buckets; c decrements through
		 * the bucket, and ends up pointing to the first element of
		 * the bucket.
		 */
		for (i = nmemb; --i >= 0;) {
			--p;
			l2[--c[tr[(*p)[indx]]]] = *p;
		}

		/* Copy the partitioned elements back to user stack */
		bcopy(l2, bot, nmemb * sizeof(u_char *));

		++indx;
		/*
		 * Sort buckets as necessary; don't sort c[0], it's the
		 * EOS character bucket, and nothing can follow EOS.
		 */
		for (i = max; i; --i) {
			if ((nmemb = c[i + 1] - (t1 = c[i])) < 2)
				continue;
			p = bot + t1;
			if (nmemb > __rspartition)
				STACKPUSH
			else
				shellsort(p, indx, nmemb, tr);
		}
		for (i = max + 1; i < NBUCKETS; ++i) {
			if ((nmemb = c[i + 1] - (t1 = c[i])) < 2)
				continue;
			p = bot + t1;
			if (nmemb > __rspartition)
				STACKPUSH
			else
				shellsort(p, indx, nmemb, tr);
		}
		/* Break out when stack is empty */
		STACKPOP
	}

	free((char *)l2);
	free((char *)stack);
	return(0);
}

/*
 * Shellsort (diminishing increment sort) from Data Structures and
 * Algorithms, Aho, Hopcraft and Ullman, 1983 Edition, page 290;
 * see also Knuth Vol. 3, page 84.  The increments are selected from
 * formula (8), page 95.  Roughly O(N^3/2).
 */
static void
shellsort(p, indx, nmemb, tr)
	register u_char **p, *tr;
	register int indx, nmemb;
{
	register u_char ch, *s1, *s2;
	register int incr, *incrp, t1, t2;

	for (incrp = __rsshell_increments; incr = *incrp++;)
		for (t1 = incr; t1 < nmemb; ++t1)
			for (t2 = t1 - incr; t2 >= 0;) {
				s1 = p[t2] + indx;
				s2 = p[t2 + incr] + indx;
				while ((ch = tr[*s1++]) == tr[*s2] && ch)
					++s2;
				if (ch > tr[*s2]) {
					s1 = p[t2];
					p[t2] = p[t2 + incr];
					p[t2 + incr] = s1;
					t2 -= incr;
				} else
					break;
			}
}
Commit	Line	Data
15637ed4 RG	1	/*-
	2	* Copyright (c) 1990 The Regents of the University of California.
	3	* All rights reserved.
	4	*
	5	* Redistribution and use in source and binary forms, with or without
	6	* modification, are permitted provided that the following conditions
	7	* are met:
	8	* 1. Redistributions of source code must retain the above copyright
	9	* notice, this list of conditions and the following disclaimer.
	10	* 2. Redistributions in binary form must reproduce the above copyright
	11	* notice, this list of conditions and the following disclaimer in the
	12	* documentation and/or other materials provided with the distribution.
	13	* 3. All advertising materials mentioning features or use of this software
	14	* must display the following acknowledgement:
	15	* This product includes software developed by the University of
	16	* California, Berkeley and its contributors.
	17	* 4. Neither the name of the University nor the names of its contributors
	18	* may be used to endorse or promote products derived from this software
	19	* without specific prior written permission.
	20	*
	21	* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
	22	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
	23	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
	24	* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
	25	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
	26	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
	27	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
	28	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
	29	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
	30	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
	31	* SUCH DAMAGE.
	32	*/
	33
	34	#if defined(LIBC_SCCS) && !defined(lint)
	35	static char sccsid[] = "@(#)radixsort.c 5.7 (Berkeley) 2/23/91";
	36	#endif /* LIBC_SCCS and not lint */
	37
	38	#include <sys/types.h>
	39	#include <limits.h>
	40	#include <stdlib.h>
	41	#include <stddef.h>
	42	#include <string.h>
	43
	44	/*
	45	* __rspartition is the cutoff point for a further partitioning instead
	46	* of a shellsort. If it changes check __rsshell_increments. Both of
	47	* these are exported, as the best values are data dependent.
	48	*/
	49	#define NPARTITION 40
	50	int __rspartition = NPARTITION;
	51	int __rsshell_increments[] = { 4, 1, 0, 0, 0, 0, 0, 0 };
	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	* all the buckets but one contain (npartitions + 1) elements and the bucket
	59	* pushed on the stack last contains the rest of the elements. In this case,
	60	* stack growth is bounded by:
	61	*
	62	* limit = (nelements / (npartitions + 1)) - 1;
	63	*
	64	* This is a very large number, 52,377,648 for the maximum 32-bit signed int.
65	*
66	* By forcing the largest bucket to be pushed on the stack first, the worst
67	* case is when all but two buckets each contain (npartitions + 1) elements,
68	* with the remaining elements split equally between the first and last
69	* buckets pushed on the stack. In this case, stack growth is bounded when:
70	*
71	* for (partition_cnt = 0; nelements > npartitions; ++partition_cnt)
72	* nelements =
73	* (nelements - (npartitions + 1) * (nbuckets - 2)) / 2;
74	* The bound is:
75	*
76	* limit = partition_cnt * (nbuckets - 1);
77	*
78	* This is a much smaller number, 4590 for the maximum 32-bit signed int.
79	*/
80	#define NBUCKETS (UCHAR_MAX + 1)
81
82	typedef struct _stack {
83	const u_char **bot;
84	int indx, nmemb;
85	} CONTEXT;
86
87	#define STACKPUSH { \
88	stackp->bot = p; \
89	stackp->nmemb = nmemb; \
90	stackp->indx = indx; \
91	++stackp; \
92	}
93	#define STACKPOP { \
94	if (stackp == stack) \
95	break; \
96	--stackp; \
97	bot = stackp->bot; \
98	nmemb = stackp->nmemb; \
99	indx = stackp->indx; \
100	}
101
102	/*
103	* A variant of MSD radix sorting; see Knuth Vol. 3, page 177, and 5.2.5,
104	* Ex. 10 and 12. Also, "Three Partition Refinement Algorithms, Paige
105	* and Tarjan, SIAM J. Comput. Vol. 16, No. 6, December 1987.
106	*
107	* This uses a simple sort as soon as a bucket crosses a cutoff point,
108	* rather than sorting the entire list after partitioning is finished.
109	* This should be an advantage.
110	*
111	* This is pure MSD instead of LSD of some number of MSD, switching to
112	* the simple sort as soon as possible. Takes linear time relative to
113	* the number of bytes in the strings.
114	*/
115	int
116	#if __STDC__
117	radixsort(const u_char *l1, int nmemb, const u_char tab, u_char endbyte)
118	#else
119	radixsort(l1, nmemb, tab, endbyte)
120	const u_char **l1;
121	register int nmemb;
122	const u_char *tab;
123	u_char endbyte;
124	#endif
125	{
126	register int i, indx, t1, t2;
127	register const u_char **l2;
128	register const u_char **p;
129	register const u_char **bot;
130	register const u_char *tr;
131	CONTEXT stack, stackp;
132	int c[NBUCKETS + 1], max;
133	u_char ltab[NBUCKETS];
134	static void shellsort();
135
136	if (nmemb <= 1)
137	return(0);
138
139	/*
140	* T1 is the constant part of the equation, the number of elements
141	* represented on the stack between the top and bottom entries.
142	* It doesn't get rounded as the divide by 2 rounds down (correct
143	* for a value being subtracted). T2, the nelem value, has to be
144	* rounded up before each divide because we want an upper bound;
145	* this could overflow if nmemb is the maximum int.
146	*/
147	t1 = ((__rspartition + 1) * (NBUCKETS - 2)) >> 1;
148	for (i = 0, t2 = nmemb; t2 > __rspartition; i += NBUCKETS - 1)
149	t2 = ((t2 + 1) >> 1) - t1;
150	if (i) {
151	if (!(stack = stackp = (CONTEXT )malloc(i sizeof(CONTEXT))))
152	return(-1);
153	} else
154	stack = stackp = NULL;
155
156	/*
157	* There are two arrays, one provided by the user (l1), and the
158	* temporary one (l2). The data is sorted to the temporary stack,
159	* and then copied back. The speedup of using index to determine
160	* which stack the data is on and simply swapping stacks back and
161	* forth, thus avoiding the copy every iteration, turns out to not
162	* be any faster than the current implementation.
163	*/
164	if (!(l2 = (const u_char *)malloc(sizeof(u_char ) * nmemb)))
165	return(-1);
166
167	/*
168	* Tr references a table of sort weights; multiple entries may
169	* map to the same weight; EOS char must have the lowest weight.
170	*/
171	if (tab)
172	tr = tab;
173	else {
174	for (t1 = 0, t2 = endbyte; t1 < t2; ++t1)
175	ltab[t1] = t1 + 1;
176	ltab[t2] = 0;
177	for (t1 = endbyte + 1; t1 < NBUCKETS; ++t1)
178	ltab[t1] = t1;
179	tr = ltab;
180	}
181
182	/* First sort is entire stack */
183	bot = l1;
184	indx = 0;
185
186	for (;;) {
187	/* Clear bucket count array */
188	bzero((char *)c, sizeof(c));
189
190	/*
191	* Compute number of items that sort to the same bucket
192	* for this index.
193	*/
194	for (p = bot, i = nmemb; --i >= 0;)
195	++c[tr[(*p++)[indx]]];
196
197	/*
198	* Sum the number of characters into c, dividing the temp
199	* stack into the right number of buckets for this bucket,
200	* this index. C contains the cumulative total of keys
201	* before and included in this bucket, and will later be
202	* used as an index to the bucket. c[NBUCKETS] contains
203	* the total number of elements, for determining how many
204	* elements the last bucket contains. At the same time
205	* find the largest bucket so it gets pushed first.
206	*/
207	for (i = max = t1 = 0, t2 = __rspartition; i <= NBUCKETS; ++i) {
208	if (c[i] > t2) {
209	t2 = c[i];
210	max = i;
211	}
212	t1 = c[i] += t1;
213	}
214
215	/*
216	* Partition the elements into buckets; c decrements through
217	* the bucket, and ends up pointing to the first element of
218	* the bucket.
219	*/
220	for (i = nmemb; --i >= 0;) {
221	--p;
222	l2[--c[tr[(p)[indx]]]] = p;
223	}
224
225	/* Copy the partitioned elements back to user stack */
226	bcopy(l2, bot, nmemb * sizeof(u_char *));
227
228	++indx;
229	/*
230	* Sort buckets as necessary; don't sort c[0], it's the
231	* EOS character bucket, and nothing can follow EOS.
232	*/
233	for (i = max; i; --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(p, indx, nmemb, tr);
241	}
242	for (i = max + 1; i < NBUCKETS; ++i) {
243	if ((nmemb = c[i + 1] - (t1 = c[i])) < 2)
244	continue;
245	p = bot + t1;
246	if (nmemb > __rspartition)
247	STACKPUSH
248	else
249	shellsort(p, indx, nmemb, tr);
250	}
251	/* Break out when stack is empty */
252	STACKPOP
253	}
254
255	free((char *)l2);
256	free((char *)stack);
257	return(0);
258	}
259
260	/*
261	* Shellsort (diminishing increment sort) from Data Structures and
262	* Algorithms, Aho, Hopcraft and Ullman, 1983 Edition, page 290;
263	* see also Knuth Vol. 3, page 84. The increments are selected from
264	* formula (8), page 95. Roughly O(N^3/2).
265	*/
266	static void
267	shellsort(p, indx, nmemb, tr)
268	register u_char *p, tr;
269	register int indx, nmemb;
270	{
271	register u_char ch, s1, s2;
272	register int incr, *incrp, t1, t2;
273
274	for (incrp = __rsshell_increments; incr = *incrp++;)
275	for (t1 = incr; t1 < nmemb; ++t1)
276	for (t2 = t1 - incr; t2 >= 0;) {
277	s1 = p[t2] + indx;
278	s2 = p[t2 + incr] + indx;
279	while ((ch = tr[s1++]) == tr[s2] && ch)
280	++s2;
281	if (ch > tr[*s2]) {
282	s1 = p[t2];
283	p[t2] = p[t2 + incr];
284	p[t2 + incr] = s1;
285	t2 -= incr;
286	} else
287	break;
288	}
289	}