home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / unix / volume25 / pty4 / part08 / radixsort.c < prev   
Encoding:
C/C++ Source or Header  |  1992-02-18  |  19.6 KB  |  546 lines
  1. /* radixsort.c, radixsort.h: linear-per-byte in-memory string sort
  2. Daniel J. Bernstein, brnstnd@nyu.edu; Keith Bostic, bostic@ucbvax.berkeley.edu.
  3. No dependencies.
  4. Requires malloc, free, bcopy, bzero. Can use -DUCHAR_MAX.
  5. 10/5/91 DJB: Made c static, added ctr.
  6. 7/23/91 DJB: Baseline. radixsort/DJB 3.0. See bottom for BSD copyright.
  7. No known patent problems.
  8.  
  9. Documentation in radixsort.3, portions based on BSD documentation.
  10.  
  11. History: I discovered adaptive distribution sort in 1987. (I haven't
  12. published a paper on it---it's too trivial for that---though I did
  13. mention it in a letter to Knuth.) It grew out of a suggestion quoted in
  14. Knuth's section on quicksort et al., volume 3, page 128: ``John McCarthy
  15. has suggested setting K \approx (u + v)/2, if all keys are known to lie
  16. between u and v.'' What's the biggest problem with MSD radix sort? You
  17. can waste lots of time considering ranges of keys that don't even exist!
  18. In distribution sort, however, you usually start by finding the minimum
  19. and maximum keys. Here McCarthy's suggestion pops in. To sort a pile of
  20. floating-point numbers, find their minimum and maximum, divide the
  21. interval evenly into n subintervals, and split the numbers into piles by
  22. interval just as in MSD radix sort. Just as in quicksort, stop splitting
  23. when a pile can be sorted efficiently by a small-scale method. That's
  24. adaptive distribution sort, and it turns out to be hellishly fast for
  25. sorting floating-point data. The credit really belongs with McCarthy---I
  26. only generalized from 2 to n.
  27.  
  28. Adaptive distribution sort can be applied to strings, too: find the
  29. first character where two strings in the pile differ, and distribute on
  30. that character. There's no fine line between adaptive distribution sort
  31. and MSD radix sort, but in any case you get a big speed boost from
  32. sorting small piles by a small-scale method. See especially Knuth,
  33. exercise 5.2.5-10.
  34.  
  35. Computer scientists should note that this method is linear in the number
  36. of bytes being sorted. Sometime in 1989, as I recall, I saw a notice
  37. that someone had discovered an o(n log n) method of sorting n integers.
  38. The method depended on all sorts of weid bit manipulations and was
  39. utterly impractical. As the integers had to be short anyway, MSD radix
  40. sort worked in O(n) time. My guess is that most computer scientists
  41. don't learn about MSD radix sort (and hence don't know that sorting is
  42. linear-per-byte) because it's widely seen as having too big a constant
  43. factor to be practical. This radixsort() is a constructive proof that
  44. the opposite is true.
  45.  
  46. I ended up sending this code to Berkeley. Keith Bostic cleaned it up,
  47. fixed a few bugs, added a shell sort for the small case, and did several
  48. helpful optimizations. radixsort() will be part of BSD 4.4. I took back
  49. his version and modified it to what you see below. Among other things, I
  50. ported it from ANSI C back to the rest of the world, cleaned up some of
  51. the comments, added a proof that part of the method actually works,
  52. added the radixsort5(), radixsort7(), and radixsort3() variants,
  53. restored the original indentation, fixed an overly conservative estimate
  54. of the necessary stack size, and plugged a memory leak.
  55. */
  56.  
  57. #define blob unsigned char /* technically, shouldn't be typedefed */
  58.  
  59. /* KB says:
  60.  * __rspartition is the cutoff point for a further partitioning instead
  61.  * of a shellsort.  If it changes check __rsshell_increments.  Both of
  62.  * these are exported, as the best values are data dependent.
  63.  */
  64. #ifndef NPARTITION
  65. #define    NPARTITION 40
  66. #endif
  67. int __rspartition = NPARTITION;
  68. int __rsshell_increments[] = { 4, 1, 0, 0, 0, 0, 0, 0 };
  69.  
  70. /* KB says:
  71.  * Shellsort (diminishing increment sort) from Data Structures and
  72.  * Algorithms, Aho, Hopcraft and Ullman, 1983 Edition, page 290;
  73.  * see also Knuth Vol. 3, page 84. The increments are selected from
  74.  * formula (8), page 95. Roughly O(N^3/2).
  75.  */
  76. static void shellsort(p,index,n,tr)
  77. register blob **p;
  78. register blob *tr;
  79. register int index;
  80. register int n;
  81. {
  82.  register blob ch;
  83.  register blob *s1;
  84.  register blob *s2;
  85.  register int incr;
  86.  register int *incrp;
  87.  register int t1;
  88.  register int t2;
  89.  
  90.  incrp = __rsshell_increments;
  91.  while (incr = *incrp++)
  92.    for (t1 = incr;t1 < n;++t1)
  93.      for (t2 = t1 - incr;t2 >= 0;)
  94.       {
  95.        s1 = p[t2] + index;
  96.        s2 = p[t2 + incr] + index;
  97.        while ((ch = tr[*s1++]) == tr[*s2] && ch)
  98.          ++s2;
  99.        if (ch > tr[*s2])
  100.     {
  101.          s1 = p[t2];
  102.          p[t2] = p[t2 + incr];
  103.          p[t2 + incr] = s1;
  104.          t2 -= incr;
  105.     }
  106.        else
  107.          break;
  108.       }
  109. }
  110.  
  111. /* KB says:
  112.  * Stackp points to context structures, where each structure schedules a
  113.  * partitioning.  Radixsort exits when the stack is empty.
  114.  *
  115.  * If the buckets are placed on the stack randomly, the worst case is when
  116.  * all the buckets but one contain (npartitions + 1) elements and the bucket
  117.  * pushed on the stack last contains the rest of the elements.  In this case,
  118.  * stack growth is bounded by:
  119.  *
  120.  *    limit = (nelements / (npartitions + 1)) - 1;
  121.  *
  122.  * This is a very large number, 52,377,648 for the maximum 32-bit signed int.
  123.  *
  124.  * By forcing the largest bucket to be pushed on the stack first, the worst
  125.  * case is when all but two buckets each contain (npartitions + 1) elements,
  126.  * with the remaining elements split equally between the first and last
  127.  * buckets pushed on the stack.  In this case, stack growth is bounded when:
  128.  *
  129.  *    for (partition_cnt = 0; nelements > npartitions; ++partition_cnt)
  130.  *        nelements =
  131.  *            (nelements - (npartitions + 1) * (nbuckets - 2)) / 2;
  132.  *
  133.  * The bound is:
  134.  *
  135.  *    limit = partition_cnt * (nbuckets - 1);
  136.  *
  137.  * This is a much smaller number, 4590 for the maximum 32-bit signed int.
  138.  
  139. Note inserted by DJB: About time I proved this... Fix the number of
  140. buckets, b. Any given pile of n elements is split into m stack piles and
  141. b - m small-scale piles which immediately disappear. We ignore the case
  142. where the pile is split into only one pile of its original size---any
  143. pile will be split into smaller piles eventually. Say the stack is left
  144. with piles of sizes p_1 ... p_m, each at least P + 1, none equal to n,
  145. while x elements, for some x from 0 to m'P, are in small-scale piles.
  146. (Here P is the cutoff.) We must have p_1 + ... + p_m + x = n. Depending
  147. on the other (p,x) constraints chosen, we define s(n) as the maximum
  148. stack size for n elements. As the subpile distributions are independent,
  149. clearly
  150.  
  151.   s(n) = max_m max_(p,x) max {s(p_1),s(p_2) + 1,...,s(p_m) + m - 1}
  152.  
  153. over all valid m and (p,x). In particular, if m > 0 then we must have
  154. n > p_1 >= P + 1, so if n <= P + 1 then the maximum is (vacuously) 0. So
  155. s(n) is monotone increasing for n <= P + 1. An easy induction shows that
  156. s(n) is in fact 0 for all n < 2P + 2.
  157.  
  158. Clearly s(n) is monotone: for n >= P + 2 we choose m = 1, p_1 = n - 1,
  159. and x = 0, and we see s(n) >= s(p_1) = s(n - 1). For m = 0, the maximum
  160. always equals 0, and for m = 1, the maximum always equals a previous
  161. s(p_1), so we have
  162.  
  163.   s(n) = max { max_{k<n} s(k) , max_{m>=2} max_(p,x) max s(p_j) + j - 1 }.
  164.  
  165. For convenience we define t(n) = max_{m>=2} max_(p,x) max s(p_j) + j - 1.
  166.  
  167. Fix n. Fix m >= 2. Consider a (p,x) achieving the maximum value of
  168. max s(p_j) + j - 1. Since p_1 >= P + 1, we have p_2 <= n - P - 1. If x
  169. isn't 0, we can move an element from one of the small-scale piles to
  170. stack pile p_2, under either of the constraints in question. This
  171. increases p_2---hence does not decrease s(p_2)---without affecting the
  172. other s(p_j). Hence there is a (p,x) with smaller x also achieving the
  173. maximum.
  174.  
  175. So consider a (p,0) achieving the maximum, and say max s(p_j) + j - 1 is
  176. achieved at j = i. If we exchange p_i and p_j while meeting the
  177. constraints, we must not be at a higher maximum; in particular,
  178. s(p_i) + j - 1 <= s(p_i) + i - 1, so j <= i.
  179.  
  180. Restrict attention the the case without constraints, NC. The choice of j
  181. is unrestricted, so in particular m <= i and hence i = m. Thus t(n)
  182. equals max_{m>=2} s(p_m) + m - 1. Since all p_j are at least P + 1, we
  183. have
  184.  
  185.   p_m + (P + 1)(m - 1) <= p_m + ... + p_1 = n
  186.   p_m <= n - (P + 1)(m - 1)
  187.   s(p_m) <= s(n - (P + 1)(m - 1))
  188.   t(n) <= max_{m>=2} s(n - (P + 1)(m - 1)) + m - 1    (NC),
  189.  
  190. and it is easy to see that for n >= (P + 1)m, we can choose p_m as
  191. n - (P + 1)(m - 1) >= P + 1 and all other p_j = P + 1, so that the bound
  192. is achieved:
  193.  
  194.   t(n) = max_{m>=2} s(n - (P + 1)(m - 1)) + m - 1  for  n/(P + 1) >= m.  (NC)
  195.  
  196. For 2 <= n/(P + 1) < b, we can choose m anywhere between 2 and
  197. f = floor(n/(P + 1)) inclusive. Now
  198.  
  199.   s(n) = max { max s(k), max_{2<=m<f} s(n - (P + 1)(m - 1)) + m - 1 } (NC)
  200.  
  201. We claim inductively that s(n) = f - 1 for any n with floor(n/(P + 1)) =
  202. f. This is true for f = 1. For larger f, s(n - (P + 1)(m - 1)) =
  203. floor((n - (P + 1)(m - 1))/(P + 1)) - 1 = f - (m - 1) - 1 = f - m, so
  204. that
  205.  
  206.   s(n) = max { max s(k), f - 1 }   (NC)
  207.  
  208. If n = (P + 1)f, then by inductive hypothesis s(k) <= f - 2, so that
  209. s(n) = f - 1 as desired. Otherwise, by (sub)induction on n, s(k) is
  210. still bounded by f - 1, so that s(n) = f - 1 always. Therefore:
  211.  
  212.   s(n) = floor(n/(P + 1)) - 1,  n >= P + 1.     (NC)
  213.  
  214. Now consider the push-largest-pile-first constraint, FC. This requires
  215. that p_1 >= p_j for all j. Hence we cannot swap p_1 with p_i. However,
  216. if i is not 1 then we can swap p_j with p_i for all j > 1, hence j <= i
  217. for all j between 2 and m, hence i is m. Thus the maximum is achieved
  218. either at i = 1 or at i = m, and we have
  219.  
  220.   t(n) = max_{m>=2} max_(p,0) max { s(p_m) + m - 1, s(p_1) }.   (FC)
  221.  
  222. Take a p vector achieving the maximum of max { s(p_m) + m - 1, s(p_1) },
  223. with m fixed. Suppose some p_j for j between 2 and m - 1 inclusive is
  224. larger than P + 1. (This is vacuous for m = 2.) Then we can subtract one
  225. from p_j and add one to p_1, preserving the constraint and not
  226. decreasing the maximum. Hence the maximum is achieved somewhere with all
  227. p_j = P + 1 for 2 <= j < m, i.e.,
  228.  
  229.   p_1 + p_m + (P + 1)(m - 2) = n.   (FC)
  230.  
  231. Furthermore, p_1 >= p_m. Hence 2p_1 >= n - (P + 1)(m - 2), and p_1 can
  232. range from ceiling((n - (P + 1)(m - 2))/2) up to n - (P + 1). Similarly,
  233. 2p_m <= n - (P + 1)(m - 2), so p_m can range from P + 1 up to
  234. floor((n - (P + 1)(m - 2))/2). The global maximum of these quantities
  235. simultaneously equals the global maximum of them individually, so if all
  236. bounds can be achieved then
  237.  
  238.   t(n) = max_{m>=2} max { s(n - (P + 1)),
  239.               s(floor((n - (P + 1)(m - 2))/2)) + m - 1 }.  (FC)
  240.  
  241. This can be achieved if n - (P + 1) >= ceiling((n - (P + 1)(m - 2))/2)
  242. and, equivalently, P + 1 <= floor((n - (P + 1)(m - 2))/2), since in that
  243. case we can choose both p_1 and p_m as stated. These reduce after some
  244. simple manipulation to n >= (P + 1)m, i.e., m <= floor(n/(P + 1)). For
  245. other m it is not possible to choose any valid p_i (exercise).
  246.  
  247. We'd like to show inductively that for all n >= 2P + 2 we have
  248.  
  249.   s(n) = u(n) with
  250.   u(n) = max_{m>=2} s(floor((n - (P + 1)(m - 2))/2)) + m - 1.   (FC,*)
  251.  
  252. To do this we need only show that the other terms of the maximum do not
  253. ``get in the way,'' i.e., that u(n) >= s(n - (P + 1)) and u(n) >= s(k)
  254. for k < n. The second half is easy: s(k) = u(k), which is at most u(n)
  255. by monotonicity of s. The first half also follows from the induction:
  256.  
  257.   s(n) = max_m s(floor((n - (P + 1)(m - 2))/2)) + m - 1
  258.   s(n - (P + 1)) = u(n - (P + 1))
  259.    = max_{m>=2} s(floor((n - (P + 1) - (P + 1)(m - 2))/2)) + m - 1
  260.    <= max_{m>=2} s(floor((n - (P + 1)(m - 2))/2)) + m - 1
  261.    = u(n)                    (FC)
  262.  
  263. again by the monotonicity of s. This proves (FC,*). For small n, i.e.,
  264. floor(n/(P + 1)) = f with 2 <= f <= b, we can choose m = f, so s(n) is
  265. at least s(floor((n - (P + 1)(f - 2))/2)) + f - 1. The floor term is at
  266. most P + 1, so s(n) >= f - 1. Furthermore, s in the constrained case is
  267. at most s in the unconstrained case, so s(n) = f - 1. For larger n, by
  268. similar logic, the maximum is attained at m = b, so we finally have
  269.  
  270.   s(n) = floor(n/(P + 1)) - 1  for n < (b + 1)(P + 1)
  271.   s(n) = s(floor((n - (P + 1)(b - 2))/2)) + b - 1  otherwise    (FC)
  272.  
  273. As in the first case s(n) is always bounded by b - 1, we can calculate
  274. a bound on s(n) by repeatedly setting n to floor(n - (P + 1)(b - 2))/2)
  275. until it is under 2P + 2 (so that s(n) = 0), counting the number of
  276. iterations necessary, and multiplying by b - 1. And that's what we
  277. wanted to prove.
  278. */
  279.  
  280. #ifndef UCHAR_MAX /* XXX: we aren't even giving a chance for a definition! */
  281. #define UCHAR_MAX 256
  282. #endif
  283. #define    NBUCKETS (UCHAR_MAX + 1)
  284.  
  285. typedef struct
  286.  {
  287.   blob **bot;
  288.   int index;
  289.   int n;
  290.  }
  291. context;
  292.  
  293. #define    STACKPUSH { \
  294.   stackp->bot = p; \
  295.   stackp->n = n; \
  296.   stackp->index = index; \
  297.   ++stackp; \
  298. }
  299.  
  300. #define    STACKPOP { \
  301.   if (stackp == stack) \
  302.     break; \
  303.   --stackp; \
  304.   bot = stackp->bot; \
  305.   n = stackp->n; \
  306.   index = stackp->index; \
  307. }
  308.  
  309. /* KB says:
  310.  * This uses a simple sort as soon as a bucket crosses a cutoff point,
  311.  * rather than sorting the entire list after partitioning is finished.
  312.  * This should be an advantage.
  313.  
  314. Note from DJB: The original comment read ``There is no strong evidence
  315. that this is an advantage.'' Depressing. Here's what I wrote in response:
  316.  
  317.    Of course it's an advantage: it has to be, I coded it that way. :-)
  318.    Seriously, I just coded the sort that way since I was following Knuth's
  319.    description of MSD to the hilt. As you can imagine, though, doing the
  320.    sort this way saves just a bit of paging of the index array. It also
  321.    means that the simple sort doesn't have to worry about crossing past
  322.    already-determined boundaries---for an average 2x gain. Trust me, it's
  323.    an advantage.
  324. */
  325.  
  326. int radixsort7(l1,n,endchar,tab,indexstart,ralloc,rfree)
  327. blob **l1;
  328. register int n;
  329. unsigned int endchar; /* could use blob, but chars are unsafe with prototypes */
  330. blob *tab;
  331. int indexstart;
  332. char *(*ralloc)();
  333. void (*rfree)();
  334. {
  335.  register int i;
  336.  register int index;
  337.  register int t1;
  338.  register int t2;
  339.  register blob **l2;
  340.  register blob **p;
  341.  register blob **bot;
  342.  register blob *tr;
  343.  context *stack;
  344.  context *stackp;
  345.  static int c[NBUCKETS + 1];
  346.  static int *ctr[NBUCKETS + 1];
  347.  int max;
  348.  blob ltab[NBUCKETS]; /* local (default) table */
  349.  
  350.  if (n <= 1)
  351.    return 0;
  352.  
  353.  /* Allocate the stack. */
  354.  t1 = (__rspartition + 1) * (NBUCKETS - 2);
  355.  for (i = 0,t2 = n;t2 > __rspartition;i += NBUCKETS - 1)
  356.    t2 = (t2 - t1)/2; /* could go negative! but that's okay */
  357.  if (!i)
  358.    stack = stackp = 0;
  359.  else
  360.    if (!(stack = stackp = (context *)ralloc(i * sizeof(context))))
  361.      return -1;
  362.  
  363.  /* KB says:
  364.   * There are two arrays, one provided by the user (l1), and the
  365.   * temporary one (l2). The data is sorted to the temporary stack,
  366.   * and then copied back. The speedup of using index to determine
  367.   * which stack the data is on and simply swapping stacks back and
  368.   * forth, thus avoiding the copy every iteration, turns out to not
  369.   * be any faster than the current implementation.
  370.   */
  371.  if (!(l2 = (blob **)ralloc(n * sizeof(blob *))))
  372.   {
  373.    rfree((char *)stackp);
  374.    return -1;
  375.   }
  376.  
  377.  /* KB says:
  378.   * tr references a table of sort weights; multiple entries may
  379.   * map to the same weight; EOS char must have the lowest weight.
  380.   */
  381.  if (tab)
  382.    tr = tab;
  383.  else
  384.   {
  385.    t2 = endchar;
  386.    for (t1 = 0;t1 < t2;++t1)
  387.      ltab[t1] = t1 + 1;
  388.    ltab[t2] = 0;
  389.    for (t1 = t2 + 1;t1 < NBUCKETS;++t1)
  390.      ltab[t1] = t1;
  391.    tr = ltab;
  392.   }
  393.  
  394.  for (t1 = 0;t1 < NBUCKETS;++t1)
  395.    ctr[t1] = c + tr[t1];
  396.  
  397.  /* First sort is entire pile. */
  398.  bot = l1;
  399.  index = indexstart;
  400.  
  401.  for (;;)
  402.   {
  403.    /* Clear the bucket count array. XXX: This isn't portable to */
  404.    /* machines where the byte representation of int 0 isn't all */
  405.    /* zeros. :-) */
  406.    bzero((char *)c,sizeof(c));
  407.  
  408.    /* Compute number of items that sort to the same bucket for this index. */
  409.    p = bot;
  410.    i = n;
  411.    while (--i >= 0)
  412.      ++*ctr[(*p++)[index]];
  413.  
  414.    /* KB says:
  415.     * Sum the number of characters into c, dividing the temp stack
  416.     * into the right number of buckets for this bucket, this index.
  417.     * c contains the cumulative total of keys before and included in
  418.     * this bucket, and will later be used as an index to the bucket. 
  419.     * c[NBUCKETS] contains the total number of elements, for determining
  420.     * how many elements the last bucket contains. At the same time, we
  421.     * find the largest bucket so it gets pushed first.
  422.     */
  423.    t2 = __rspartition;
  424.    max = t1 = 0;
  425.    for (i = 0;i <= NBUCKETS;++i)
  426.     {
  427.      if (c[i] > t2)
  428.       {
  429.        t2 = c[i];
  430.        max = i;
  431.       }
  432.      t1 = c[i] += t1;
  433.     }
  434.  
  435.    /* Partition the elements into buckets; c decrements through the */
  436.    /* bucket, and ends up pointing to the first element of the bucket. */
  437.    i = n;
  438.    while (--i >= 0)
  439.     {
  440.      --p;
  441.      l2[--*ctr[(*p)[index]]] = *p;
  442.     }
  443.  
  444.    /* Copy the partitioned elements back to the user stack. */
  445.    bcopy(l2,bot,n * sizeof(blob *));
  446.  
  447.    ++index;
  448.    /* KB says:
  449.     * Sort buckets as necessary; don't sort c[0], it's the
  450.     * EOS character bucket, and nothing can follow EOS.
  451.     */
  452.    for (i = max;i;--i)
  453.     {
  454.      if ((n = c[i + 1] - (t1 = c[i])) < 2)
  455.        continue;
  456.      p = bot + t1;
  457.      if (n > __rspartition)
  458.        STACKPUSH
  459.      else
  460.        shellsort(p,index,n,tr);
  461.     }
  462.    for (i = max + 1;i < NBUCKETS;++i)
  463.     {
  464.      if ((n = c[i + 1] - (t1 = c[i])) < 2)
  465.        continue;
  466.      p = bot + t1;
  467.      if (n > __rspartition)
  468.        STACKPUSH
  469.      else
  470.        shellsort(p,index,n,tr);
  471.     }
  472.    /* Break out when stack is empty */
  473.    STACKPOP
  474.   }
  475.  
  476.  rfree((char *)l2);
  477.  rfree((char *)stack);
  478.  return 0;
  479. }
  480.  
  481. extern char *malloc();
  482. extern void free();
  483.  
  484. int radixsort5(l1,n,endchar,tab,indexstart)
  485. blob **l1; register int n; unsigned int endchar; blob *tab; int indexstart;
  486. {
  487.  return radixsort7(l1,n,endchar,tab,indexstart,malloc,free);
  488. }
  489.  
  490. int radixsort4(l1,n,endchar,tab)
  491. blob **l1; register int n; unsigned int endchar; blob *tab;
  492. {
  493.  return radixsort7(l1,n,endchar,tab,0,malloc,free);
  494. }
  495.  
  496. int radixsort(l1,n,tab,endchar)
  497. blob **l1; register int n; blob *tab; unsigned int endchar;
  498. {
  499.  return radixsort7(l1,n,endchar,tab,0,malloc,free);
  500. }
  501.  
  502. int radixsort3(l1,n,endchar)
  503. blob **l1; register int n; unsigned int endchar; 
  504. {
  505.  return radixsort7(l1,n,endchar,(blob *) 0,0,malloc,free);
  506. }
  507.  
  508. /*- BSD copyright says:
  509.  * Copyright (c) 1990 The Regents of the University of California.
  510.  * All rights reserved.
  511.  *
  512.  * Redistribution and use in source and binary forms, with or without
  513.  * modification, are permitted provided that the following conditions
  514.  * are met:
  515.  * 1. Redistributions of source code must retain the above copyright
  516.  *    notice, this list of conditions and the following disclaimer.
  517.  * 2. Redistributions in binary form must reproduce the above copyright
  518.  *    notice, this list of conditions and the following disclaimer in the
  519.  *    documentation and/or other materials provided with the distribution.
  520.  * 3. All advertising materials mentioning features or use of this software
  521.  *    must display the following acknowledgement:
  522.  *    This product includes software developed by the University of
  523.  *    California, Berkeley and its contributors.
  524.  * 4. Neither the name of the University nor the names of its contributors
  525.  *    may be used to endorse or promote products derived from this software
  526.  *    without specific prior written permission.
  527.  *
  528.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  529.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  530.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  531.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  532.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  533.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  534.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  535.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  536.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  537.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  538.  * SUCH DAMAGE.
  539.  */
  540.  
  541. /* BSD sccs says:
  542. static char sccsid[] = "@(#)radixsort.c  5.7 (Berkeley) 2/23/91";
  543.  
  544. but this is (heavily) modified.
  545. */
  546.