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