added radix sort

diff --git a/test/ksort_test.cc b/test/ksort_test.cc
index 2cf8b8a4563c1d862603f3bc555ca71a59622ab2..836a52eef46b847e8f1f2dda9493547ed2692515 100644 (file)
--- a/test/ksort_test.cc
+++ b/test/ksort_test.cc
@@ -568,6 +568,92 @@ int * tmpArray;
  * END OF PAUL'S IMPLEMENTATION *
  ********************************/
 
+#define rstype_t uint32_t
+#define rskey(x) (x)
+
+#define RS_MIN_SIZE 63
+
+typedef struct {
+       rstype_t *b, *e;
+} bucket_t;
+
+static inline void rs_insertsort(rstype_t *s, rstype_t *t)
+{
+       rstype_t *i, *j, swap_tmp;
+       for (i = s + 1; i < t; ++i)
+               for (j = i; j > s && rskey(*j) < rskey(*(j-1)); --j) {
+                       swap_tmp = *j; *j = *(j-1); *(j-1) = swap_tmp;
+               }
+}
+
+void rs_combsort(size_t n, rstype_t a[])
+{
+       const double shrink_factor = 1. / 1.2473309501039786540366528676643;
+       int do_swap;
+       size_t gap = n;
+       rstype_t tmp, *i, *j;
+       do {
+               if (gap > 2) {
+                       gap = (size_t)(gap * shrink_factor);
+                       if (gap == 9 || gap == 10) gap = 11;
+               }
+               do_swap = 0;
+               for (i = a; i < a + n - gap; ++i) {
+                       j = i + gap;
+                       if (rskey(*j) < rskey(*i)) {
+                               tmp = *i; *i = *j; *j = tmp;
+                               do_swap = 1;
+                       }
+               }
+       } while (do_swap || gap > 2);
+       if (gap != 1) rs_insertsort(a, a + n);
+}
+
+void rs_classify(rstype_t *beg, rstype_t *end, int n_bits, int s, bucket_t *b)
+{
+       rstype_t *i, tmp;
+       int m = (1<<n_bits) - 1;
+       bucket_t *k, *l, *be;
+
+       be = b + (1<<n_bits);
+       for (k = b; k != be; ++k) k->b = k->e = beg;
+       for (i = beg; i != end; ++i) ++b[rskey(*i)>>s&m].e;
+       if (b[0].e == end) return; // no need to sort
+       for (k = b + 1; k != be; ++k)
+               k->e += (k-1)->e - beg, k->b = (k-1)->e;
+       for (k = b; k != be;) {
+               if (k->b == k->e) { ++k; continue; }
+               l = b + (rskey(*k->b)>>s&m);
+               if (k == l) { ++k->b; continue; }
+               while (b + (rskey(*l->b)>>s&m) == l) ++l->b;
+               tmp = *l->b; *l->b++ = *k->b; *k->b = tmp;
+       }
+       for (k = b + 1; k != be; ++k) k->b = (k-1)->e;
+       b->b = beg;
+}
+
+void rs_sort(rstype_t *beg, rstype_t *end, int n_bits, int s)
+{
+       if (end - beg < 2) { // already sorted
+               return;
+       } else if (end - beg > RS_MIN_SIZE) {
+               bucket_t *b;
+               int i;
+               b = (bucket_t*)malloc(sizeof(bucket_t) * (1<<n_bits));
+               rs_classify(beg, end, n_bits, s, b);
+               if (s) {
+                       s = s > n_bits? s - n_bits : 0;
+                       for (i = 0; i != 1<<n_bits; ++i)
+                               if (b[i].e > b[i].b + 1) rs_sort(b[i].b, b[i].e, n_bits, s);
+               }
+               free(b);
+       } else rs_combsort(end - beg, beg);
+}
+
+/*************************
+ *** END OF RADIX SORT ***
+ *************************/
+
 struct intcmp_t {
        inline int operator() (int a, int b) const {
                return a < b? -1 : a > b? 1 : 0;
@@ -593,6 +679,23 @@ int main(int argc, char *argv[])
        temp = (int*)malloc(sizeof(int) * N);
        array = (int*)malloc(sizeof(int) * N);
 
+       srand48(11);
+       for (i = 0; i < N; ++i) array[i] = (int)lrand48();
+       t1 = clock();
+       rs_sort((uint32_t*)array, (uint32_t*)array + N, 8, 24);
+       t2 = clock();
+       fprintf(stderr, "radix sort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
+       for (i = 0; i < N-1; ++i) {
+               if (array[i] > array[i+1]) {
+                       fprintf(stderr, "Bug in radix sort!\n");
+                       exit(1);
+               }
+       }
+       t1 = clock();
+       rs_sort((uint32_t*)array, (uint32_t*)array + N, 8, 24);
+       t2 = clock();
+       fprintf(stderr, "radix sort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
+
        srand48(11);
        for (i = 0; i < N; ++i) array[i] = (int)lrand48();
        t1 = clock();