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