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.3 (Berkeley) 10/13/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 * all the buckets but one contain (NPARTITION + 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 * (nelements / (npartitions + 1)) - 1 63 * 64 * This is a very large number. By forcing the largest bucket to be pushed 65 * on the stack first the worst case is when all but two buckets each contain 66 * (NPARTITION + 1) elements, with the remaining elements split equally between 67 * the first and last buckets pushed on the stack. In this case, stack growth 68 * is bounded when: 69 * 70 * for (partition_cnt = 0; nelements > npartitions; ++partition_cnt) 71 * nelements = 72 * (nelements - (npartitions + 1) * (nbuckets - 2)) / 2; 73 * The bound is: 74 * 75 * limit = partition_cnt * (nbuckets - 1); 76 * 77 * This is a much smaller number. 78 */ 79 typedef struct _stack { 80 u_char **bot; 81 int indx, nmemb; 82 } CONTEXT; 83 84 #define STACKPUSH { \ 85 stackp->bot = p; \ 86 stackp->nmemb = nmemb; \ 87 stackp->indx = indx; \ 88 ++stackp; \ 89 } 90 #define STACKPOP { \ 91 if (stackp == stack) \ 92 break; \ 93 --stackp; \ 94 bot = stackp->bot; \ 95 nmemb = stackp->nmemb; \ 96 indx = stackp->indx; \ 97 } 98 99 /* 100 * A variant of MSD radix sorting; see Knuth Vol. 3, page 177, and 5.2.5, 101 * Ex. 10 and 12. Also, "Three Partition Refinement Algorithms, Paige 102 * and Tarjan, SIAM J. Comput. Vol. 16, No. 6, December 1987. 103 * 104 * This uses a simple sort as soon as a bucket crosses a cutoff point, 105 * rather than sorting the entire list after partitioning is finished. 106 * This should be an advantage. 107 * 108 * This is pure MSD instead of LSD of some number of MSD, switching to 109 * the simple sort as soon as possible. Takes linear time relative to 110 * the number 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 *stack, *stackp; 119 int c[NCHARS + 1], max; 120 u_char ltab[NCHARS]; 121 122 if (nmemb <= 1) 123 return(0); 124 125 /* 126 * T1 is the constant part of the equation, the number of elements 127 * represented on the stack between the top and bottom entries. 128 * It doesn't get rounded as the divide by 2 rounds down (correct 129 * for a value being subtracted). T2, the nelem value, has to be 130 * rounded up before each divide because we want an upper bound; 131 * this could overflow if nmemb is the maximum int. 132 */ 133 t1 = ((__rspartition + 1) * (UCHAR_MAX - 2)) >> 1; 134 for (i = 0, t2 = nmemb; t2 > __rspartition; i += UCHAR_MAX - 1) 135 t2 = (++t2 >> 1) - t1; 136 if (i) { 137 if (!(stack = stackp = (CONTEXT *)malloc(i * sizeof(CONTEXT)))) 138 return(-1); 139 } else 140 stack = stackp = NULL; 141 142 /* 143 * There are two arrays, one provided by the user (l1), and the 144 * temporary one (l2). The data is sorted to the temporary stack, 145 * and then copied back. The speedup of using index to determine 146 * which stack the data is on and simply swapping stacks back and 147 * forth, thus avoiding the copy every iteration, turns out to not 148 * be any faster than the current implementation. 149 */ 150 if (!(l2 = (u_char **)malloc(sizeof(u_char *) * nmemb))) 151 return(-1); 152 153 /* 154 * Tr references a table of sort weights; multiple entries may 155 * map to the same weight; EOS char must have the lowest weight. 156 */ 157 if (tab) 158 tr = tab; 159 else { 160 tr = ltab; 161 for (t1 = 0, t2 = endbyte; t1 < t2; ++t1) 162 tr[t1] = t1 + 1; 163 tr[t2] = 0; 164 for (t1 = endbyte + 1; t1 < NCHARS; ++t1) 165 tr[t1] = t1; 166 } 167 168 /* First sort is entire stack */ 169 bot = l1; 170 indx = 0; 171 172 for (;;) { 173 /* Clear bucket count array */ 174 bzero((char *)c, sizeof(c)); 175 176 /* 177 * Compute number of items that sort to the same bucket 178 * for this index. 179 */ 180 for (p = bot, i = nmemb; i--;) 181 ++c[tr[(*p++)[indx]]]; 182 183 /* 184 * Sum the number of characters into c, dividing the temp 185 * stack into the right number of buckets for this bucket, 186 * this index. C contains the cumulative total of keys 187 * before and included in this bucket, and will later be 188 * used as an index to the bucket. c[NCHARS] contains 189 * the total number of elements, for determining how many 190 * elements the last bucket contains. At the same time 191 * find the largest bucket so it gets handled first. 192 */ 193 for (i = 1, t2 = -1; i <= NCHARS; ++i) { 194 if ((t1 = c[i - 1]) > t2) { 195 t2 = t1; 196 max = i; 197 } 198 c[i] += t1; 199 } 200 201 /* 202 * Partition the elements into buckets; c decrements 203 * through the bucket, and ends up pointing to the 204 * first element of the bucket. 205 */ 206 for (i = nmemb; i--;) { 207 --p; 208 l2[--c[tr[(*p)[indx]]]] = *p; 209 } 210 211 /* Copy the partitioned elements back to user stack */ 212 bcopy(l2, bot, nmemb * sizeof(u_char *)); 213 214 ++indx; 215 /* 216 * Sort buckets as necessary; don't sort c[0], it's the 217 * EOS character bucket, and nothing can follow EOS. 218 */ 219 for (i = max; i; --i) { 220 if ((nmemb = c[i + 1] - (t1 = c[i])) < 2) 221 continue; 222 p = bot + t1; 223 if (nmemb > __rspartition) 224 STACKPUSH 225 else 226 SHELLSORT 227 } 228 for (i = max + 1; i < NCHARS; ++i) { 229 if ((nmemb = c[i + 1] - (t1 = c[i])) < 2) 230 continue; 231 p = bot + t1; 232 if (nmemb > __rspartition) 233 STACKPUSH 234 else 235 SHELLSORT 236 } 237 /* Break out when stack is empty */ 238 STACKPOP 239 } 240 241 free((char *)l2); 242 free((char *)stack); 243 return(0); 244 } 245