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