c_cpp/lib/klib/test/ksort_test.cc

   1 #include <stdlib.h>
   2 #include <string.h>
   3 #include <stdio.h>
   4 #include <time.h>
   5 #include <algorithm>
   6
   7 #include "ksort.h"
   8 KSORT_INIT_GENERIC(int)
   9
  10 using namespace std;
  11
  12 /**********************************
  13  * BEGIN OF PAUL'S IMPLEMENTATION *
  14  **********************************/
  15
  16 /* Attractive Chaos: I have added inline where necessary. */
  17
  18 /*
  19 Copyright (c) 2004 Paul Hsieh
  20 All rights reserved.
  21
  22 Redistribution and use in source and binary forms, with or without
  23 modification, are permitted provided that the following conditions are met:
  24
  25     Redistributions of source code must retain the above copyright notice,
  26     this list of conditions and the following disclaimer.
  27
  28     Redistributions in binary form must reproduce the above copyright notice,
  29     this list of conditions and the following disclaimer in the documentation
  30     and/or other materials provided with the distribution.
  31
  32     Neither the name of sorttest nor the names of its contributors may be
  33     used to endorse or promote products derived from this software without
  34     specific prior written permission.
  35
  36 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  37 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  38 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  39 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  40 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  41 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  42 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  43 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  44 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  45 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  46 POSSIBILITY OF SUCH DAMAGE.
  47 */
  48
  49 /*
  50
  51 Recommended flags:
  52 ------------------
  53
  54 Intel C/C++:
  55 icl /O2 /G6 /Qaxi /Qxi /Qip sorttest.c
  56
  57 WATCOM C/C++:
  58 wcl386 /otexan /6r sorttest.c
  59
  60 GCC:
  61 gcc -O3 -mcpu=athlon-xp -march=athlon-xp sorttest.c
  62
  63 MSVC:
  64 cl /O2 /Ot /Og /G6 sorttest.c
  65
  66 */
  67
  68 static inline void sort2 (int * numbers) {
  69 int tmp;
  70
  71     if (numbers[0] <= numbers[1]) return;
  72     tmp = numbers[0];
  73     numbers[0] = numbers[1];
  74     numbers[1] = tmp;
  75 }
  76
  77 static inline void sort3 (int * numbers) {
  78 int tmp;
  79
  80     if (numbers[0] <= numbers[1]) {
  81         if (numbers[1] <= numbers[2]) return;
  82         if (numbers[2] <= numbers[0]) {
  83             tmp = numbers[0];
  84             numbers[0] = numbers[2];
  85             numbers[2] = numbers[1];
  86             numbers[1] = tmp;
  87             return;
  88         }
  89         tmp = numbers[1];
  90     } else {
  91         tmp = numbers[0];
  92         if (numbers[0] <= numbers[2]) {
  93             numbers[0] = numbers[1];
  94             numbers[1] = tmp;
  95             return;
  96         }
  97         if (numbers[2] <= numbers[1]) {
  98             numbers[0] = numbers[2];
  99             numbers[2] = tmp;
 100             return;
 101         }
 102         numbers[0] = numbers[1];
 103     }
 104     numbers[1] = numbers[2];
 105     numbers[2] = tmp;
 106 }
 107
 108 static inline void sort4 (int * num) {
 109 int tmp;
 110   if (num[0] < num[1]) {
 111     if (num[1] < num[2]) {
 112       if (num[1] < num[3]) {
 113         if (num[2] >= num[3]) {
 114           tmp = num[2];
 115           num[2] = num[3];
 116           num[3] = tmp;
 117         }
 118       } else {
 119         tmp = num[1];
 120         if (num[0] < num[3]) {
 121           num[1] = num[3];
 122         } else {
 123           num[1] = num[0];
 124           num[0] = num[3];
 125         }
 126         num[3] = num[2];
 127         num[2] = tmp;
 128       }
 129     } else {
 130       if (num[0] < num[2]) {
 131         if (num[2] < num[3]) {
 132           if (num[1] < num[3]) {
 133             tmp = num[1];
 134           } else {
 135             tmp = num[3];
 136             num[3] = num[1];
 137           }
 138           num[1] = num[2];
 139           num[2] = tmp;
 140         } else {
 141           if (num[0] < num[3]) {
 142             tmp = num[3];
 143           } else {
 144             tmp = num[0];
 145             num[0] = num[3];
 146           }
 147           num[3] = num[1];
 148           num[1] = tmp;
 149         }
 150       } else {
 151         if (num[0] < num[3]) {
 152           tmp = num[0];
 153           num[0] = num[2];
 154           if (num[1] < num[3]) {
 155             num[2] = num[1];
 156           } else {
 157             num[2] = num[3];
 158             num[3] = num[1];
 159           }
 160           num[1] = tmp;
 161         } else {
 162           if (num[2] < num[3]) {
 163             tmp = num[0];
 164             num[0] = num[2];
 165             num[2] = tmp;
 166             tmp = num[1];
 167             num[1] = num[3];
 168           } else {
 169             tmp = num[1];
 170             num[1] = num[2];
 171             num[2] = num[0];
 172             num[0] = num[3];
 173           }
 174           num[3] = tmp;
 175         }
 176       }
 177     }
 178   } else {
 179     tmp = num[0];
 180     if (tmp < num[2]) {
 181       if (tmp < num[3]) {
 182         num[0] = num[1];
 183         num[1] = tmp;
 184         if (num[2] >= num[3]) {
 185           tmp = num[2];
 186           num[2] = num[3];
 187           num[3] = tmp;
 188         }
 189       } else {
 190         if (num[1] < num[3]) {
 191           num[0] = num[1];
 192           num[1] = num[3];
 193         } else {
 194           num[0] = num[3];
 195         }
 196         num[3] = num[2];
 197         num[2] = tmp;
 198       }
 199     } else {
 200       if (num[1] < num[2]) {
 201         if (num[2] < num[3]) {
 202           num[0] = num[1];
 203           num[1] = num[2];
 204           if (tmp < num[3]) {
 205             num[2] = tmp;
 206           } else {
 207             num[2] = num[3];
 208             num[3] = tmp;
 209           }
 210         } else {
 211           if (num[1] < num[3]) {
 212             num[0] = num[1];
 213             num[1] = num[3];
 214           } else {
 215             num[0] = num[3];
 216           }
 217           num[3] = tmp;
 218         }
 219       } else {
 220         if (num[1] < num[3]) {
 221           num[0] = num[2];
 222           if (tmp < num[3]) {
 223             num[2] = tmp;
 224           } else {
 225             num[2] = num[3];
 226             num[3] = tmp;
 227           }
 228         } else {
 229           if (num[2] < num[3]) {
 230             num[0] = num[2];
 231             num[2] = num[1];
 232             num[1] = num[3];
 233             num[3] = tmp;
 234           } else {
 235             num[0] = num[3];
 236             num[3] = tmp;
 237             tmp = num[1];
 238             num[1] = num[2];
 239             num[2] = tmp;
 240           }
 241         }
 242       }
 243     }
 244   }
 245 }
 246
 247 static inline void sortAlt2 (int * numbers, int * altNumbers) {
 248     if (numbers[0] <= numbers[1]) {
 249         altNumbers[0] = numbers[0];
 250         altNumbers[1] = numbers[1];
 251     } else {
 252         altNumbers[0] = numbers[1];
 253         altNumbers[1] = numbers[0];
 254     }
 255 }
 256
 257 static inline void sortAlt3 (int * numbers, int * altNumbers) {
 258     if (numbers[0] <= numbers[1]) {
 259         if (numbers[1] <= numbers[2]) {
 260             altNumbers[0] = numbers[0];
 261             altNumbers[1] = numbers[1];
 262             altNumbers[2] = numbers[2];
 263         } else if (numbers[2] <= numbers[0]) {
 264             altNumbers[0] = numbers[2];
 265             altNumbers[1] = numbers[0];
 266             altNumbers[2] = numbers[1];
 267         } else {
 268             altNumbers[0] = numbers[0];
 269             altNumbers[1] = numbers[2];
 270             altNumbers[2] = numbers[1];
 271         }
 272     } else {
 273         if (numbers[0] <= numbers[2]) {
 274             altNumbers[0] = numbers[1];
 275             altNumbers[1] = numbers[0];
 276             altNumbers[2] = numbers[2];
 277         } else if (numbers[2] <= numbers[1]) {
 278             altNumbers[0] = numbers[2];
 279             altNumbers[1] = numbers[1];
 280             altNumbers[2] = numbers[0];
 281         } else {
 282             altNumbers[0] = numbers[1];
 283             altNumbers[1] = numbers[2];
 284             altNumbers[2] = numbers[0];
 285         }
 286     }
 287 }
 288
 289 /*
 290  *  Insert Sort
 291  */
 292
 293 inline void insertSort (int numbers[], int qty) {
 294 int i, j, idx, q4;
 295 int tmp;
 296
 297     if (qty <= 4) {
 298         if (qty == 4) sort4 (numbers);
 299         else if (qty == 3) sort3 (numbers);
 300         else if (qty == 2) sort2 (numbers);
 301         return;
 302     }
 303
 304     q4 = qty - 4;
 305
 306     for (i=0; i < q4; i++) {
 307         idx = i;
 308         for (j=i+1; j < qty; j++) {
 309             if (numbers[j] < numbers[idx]) idx = j;
 310         }
 311         if (idx != i) {
 312             tmp = numbers[idx];
 313             numbers[idx] = numbers[i];
 314             numbers[i] = tmp;
 315         }
 316     }
 317
 318     sort4 (numbers + q4);
 319 }
 320
 321 /*
 322  *  Heap Sort
 323  */
 324
 325 /* Assure the heap property for entries from top to last */
 326 static void siftDown (int numbers[], int top, int last) {
 327 int tmp = numbers[top];
 328 int maxIdx = top;
 329
 330     while (last >= (maxIdx += maxIdx)) {
 331
 332         /* This is where the comparison occurrs and where a sufficiently
 333            good compiler can use a computed conditional result rather
 334            than using control logic. */
 335         if (maxIdx != last && numbers[maxIdx] < numbers[maxIdx + 1]) maxIdx++;
 336
 337         if (tmp >= numbers[maxIdx]) break;
 338         numbers[top] = numbers[maxIdx];
 339         top = maxIdx;
 340     }
 341     numbers[top] = tmp;
 342 }
 343
 344 /* Peel off the top siftDown operation since its parameters are trivial to
 345    fill in directly (and this saves us some moves.) */
 346 static void siftDown0 (int numbers[], int last) {
 347 int tmp;
 348
 349     if (numbers[0] < numbers[1]) {
 350         tmp = numbers[1];
 351         numbers[1] = numbers[0];
 352         siftDown (numbers, 1, last);
 353     } else {
 354         tmp = numbers[0];
 355     }
 356     numbers[0] = numbers[last];
 357     numbers[last] = tmp;
 358 }
 359
 360 void heapSort (int numbers[], int qty) {
 361 int i;
 362
 363     if (qty <= 4) {
 364         if (qty == 4) sort4 (numbers);
 365         else if (qty == 3) sort3 (numbers);
 366         else if (qty == 2) sort2 (numbers);
 367         return;
 368     }
 369
 370     i = qty / 2;
 371     /* Enforce the heap property for each position in the tree */
 372     for (  qty--; i >  0; i--) siftDown  (numbers, i, qty);
 373     for (i = qty; i > 0; i--) siftDown0 (numbers, i);
 374 }
 375
 376 /*
 377  *  Quick Sort
 378  */
 379
 380 static int medianOf3 (int * numbers, int i, int j) {
 381 int tmp;
 382
 383     if (numbers[0] <= numbers[i]) {
 384         if (numbers[j] <= numbers[0]) return numbers[0]; /* j 0 i */
 385         if (numbers[i] <= numbers[j]) j = i;             /* 0 i j */
 386                                                          /* 0 j i */
 387     } else {
 388         if (numbers[0] <= numbers[j]) return numbers[0]; /* i 0 j */
 389         if (numbers[j] <= numbers[i]) j = i;             /* j i 0 */
 390                                                          /* i j 0 */
 391     }
 392     tmp = numbers[j];
 393     numbers[j] = numbers[0];
 394     numbers[0] = tmp;
 395     return tmp;
 396 }
 397
 398 static void quickSortRecurse (int * numbers, int left, int right) {
 399 int pivot, lTmp, rTmp;
 400
 401     qsrStart:;
 402
 403 #if defined(__GNUC__)
 404     if (right <= left + 8) {
 405         insertSort (numbers + left, right - left + 1);
 406         return;
 407     }
 408 #else
 409     if (right <= left + 3) {
 410         if (right == left + 1) {
 411             sort2 (numbers + left);
 412         } else if (right == left + 2) {
 413             sort3 (numbers + left);
 414         } else if (right == left + 3) {
 415             sort4 (numbers + left);
 416         }
 417         return;
 418     }
 419 #endif
 420
 421     lTmp = left;
 422     rTmp = right;
 423
 424     pivot = medianOf3 (numbers + left, (right-left) >> 1, right-1-left);
 425
 426     goto QStart;
 427     while (1) {
 428         do {
 429             right--;
 430             if (left >= right) goto QEnd;
 431             QStart:;
 432         } while (numbers[right] > pivot);
 433         numbers[left] = numbers[right];
 434         do {
 435             left++;
 436             if (left >= right) {
 437                 left = right;
 438                 goto QEnd;
 439             }
 440         } while (numbers[ left] < pivot);
 441         numbers[right] = numbers[left];
 442     }
 443     QEnd:;
 444     numbers[left] = pivot;
 445
 446     /* Only recurse the smaller partition */
 447
 448     if (left-1 - lTmp <= rTmp - left - 1) {
 449         if (lTmp < left) quickSortRecurse (numbers,   lTmp, left-1);
 450
 451         /* Set up for larger partition */
 452         left++;
 453         right = rTmp;
 454     } else {
 455         if (rTmp > left) quickSortRecurse (numbers, left+1,   rTmp);
 456
 457         /* Set up for larger partition */
 458         right = left - 1;
 459         left = lTmp;
 460     }
 461
 462     /* Rerun with larger partition (recursion not required.) */
 463     goto qsrStart;
 464 }
 465
 466 void quickSort (int numbers[], int qty) {
 467     if (qty < 2) return;
 468     quickSortRecurse (numbers, 0, qty - 1);
 469 }
 470
 471 /*
 472  *  Merge Sort
 473  */
 474
 475 static void mergesortInPlace (int * numbers, int * altNumbers, int qty);
 476
 477 /* Perform mergesort, but store results in altNumbers */
 478
 479 static void mergesortExchange (int * numbers, int * altNumbers, int qty) {
 480 int half, i0, i1, i;
 481
 482     if (qty == 2) {
 483         sortAlt2 (numbers, altNumbers);
 484         return;
 485     }
 486     if (qty == 3) {
 487         sortAlt3 (numbers, altNumbers);
 488         return;
 489     }
 490
 491     half = (qty + 1)/2;
 492
 493     mergesortInPlace (numbers, altNumbers, half);
 494     mergesortInPlace (numbers + half, altNumbers, qty - half);
 495
 496     i0 = 0; i1 = half;
 497
 498     for (i=0; i < qty; i++) {
 499         if (i1 >= qty || (i0 < half && numbers[i0] < numbers[i1])) {
 500             altNumbers[i] = numbers[i0];
 501             i0++;
 502         } else {
 503             altNumbers[i] = numbers[i1];
 504             i1++;
 505         }
 506     }
 507 }
 508
 509 /* Perform mergesort and store results in numbers */
 510
 511 static void mergesortInPlace (int * numbers, int * altNumbers, int qty) {
 512 int half, i0, i1, i;
 513
 514 #if 0
 515     if (qty == 2) {
 516         sort2 (numbers);
 517         return;
 518     }
 519     if (qty == 3) {
 520         sort3 (numbers);
 521         return;
 522     }
 523     if (qty == 4) {
 524         sort4 (numbers);
 525         return;
 526     }
 527 #else
 528     if (qty <= 12) {
 529         insertSort (numbers, qty);
 530         return;
 531     }
 532 #endif
 533
 534     half = (qty + 1)/2;
 535
 536     mergesortExchange (numbers, altNumbers, half);
 537     mergesortExchange (numbers + half, altNumbers + half, qty - half);
 538
 539     i0 = 0; i1 = half;
 540
 541     for (i=0; i < qty; i++) {
 542         if (i1 >= qty || (i0 < half && altNumbers[i0] < altNumbers[i1])) {
 543             numbers[i] = altNumbers[i0];
 544             i0++;
 545         } else {
 546             numbers[i] = altNumbers[i1];
 547             i1++;
 548         }
 549     }
 550 }
 551
 552 #include <stdlib.h>
 553
 554 void mergeSort (int numbers[], int qty) {
 555 int * tmpArray;
 556
 557     if (qty <= 12) {
 558         insertSort (numbers, qty);
 559         return;
 560     }
 561
 562     tmpArray = (int *) malloc (qty * sizeof (int));
 563     mergesortInPlace (numbers, tmpArray, qty);
 564     free (tmpArray);
 565 }
 566
 567 /********************************
 568  * END OF PAUL'S IMPLEMENTATION *
 569  ********************************/
 570
 571 /*************************************************
 572  *** Implementation 1: faster on sorted arrays ***
 573  *************************************************/
 574
 575 #define rstype_t unsigned
 576 #define rskey(x) (x)
 577
 578 #define RS_MIN_SIZE 64
 579
 580 typedef struct {
 581         rstype_t *b, *e;
 582 } rsbucket_t;
 583
 584 void rs_sort(rstype_t *beg, rstype_t *end, int n_bits, int s)
 585 {
 586         rstype_t *i;
 587         int size = 1<<n_bits, m = size - 1;
 588         rsbucket_t *k, b[size], *be = b + size;
 589
 590         for (k = b; k != be; ++k) k->b = k->e = beg;
 591         for (i = beg; i != end; ++i) ++b[rskey(*i)>>s&m].e;
 592         for (k = b + 1; k != be; ++k)
 593                 k->e += (k-1)->e - beg, k->b = (k-1)->e;
 594         for (k = b; k != be;) {
 595                 if (k->b != k->e) {
 596                         rsbucket_t *l;
 597                         if ((l = b + (rskey(*k->b)>>s&m)) != k) {
 598                                 rstype_t tmp = *k->b, swap;
 599                                 do {
 600                                         swap = tmp; tmp = *l->b; *l->b++ = swap;
 601                                         l = b + (rskey(tmp)>>s&m);
 602                                 } while (l != k);
 603                                 *k->b++ = tmp;
 604                         } else ++k->b;
 605                 } else ++k;
 606         }
 607         for (b->b = beg, k = b + 1; k != be; ++k) k->b = (k-1)->e;
 608         if (s) {
 609                 s = s > n_bits? s - n_bits : 0;
 610                 for (k = b; k != be; ++k)
 611                         if (k->e - k->b > RS_MIN_SIZE) rs_sort(k->b, k->e, n_bits, s);
 612                         else if (k->e - k->b > 1)
 613                                 for (i = k->b + 1; i < k->e; ++i)
 614                                         if (rskey(*i) < rskey(*(i - 1))) {
 615                                                 rstype_t *j, tmp = *i;
 616                                                 for (j = i; j > k->b && rskey(tmp) < rskey(*(j-1)); --j)
 617                                                         *j = *(j - 1);
 618                                                 *j = tmp;
 619                                         }
 620         }
 621 }
 622
 623 /*************************************************
 624  *** Implementation 2: faster on random arrays ***
 625  *************************************************/
 626
 627 static inline void rs_insertsort(rstype_t *s, rstype_t *t)
 628 {
 629         rstype_t *i;
 630         for (i = s + 1; i < t; ++i) {
 631                 if (rskey(*i) < rskey(*(i - 1))) {
 632                         rstype_t *j, tmp = *i;
 633                         for (j = i; j > s && rskey(tmp) < rskey(*(j-1)); --j)
 634                                 *j = *(j - 1);
 635                         *j = tmp;
 636                 }
 637         }
 638 }
 639 /*
 640 void rs_sort2(rstype_t *beg, rstype_t *end, int n_bits, int s)
 641 {
 642         int j, size = 1<<n_bits, m = size - 1;
 643         unsigned long c[size];
 644         rstype_t *i, *b[size], *e[size];
 645
 646         for (j = 0; j < size; ++j) c[j] = 0;
 647         for (i = beg; i != end; ++i) ++c[rskey(*i)>>s&m];
 648         b[0] = e[0] = beg;
 649         for (j = 1; j != size; ++j) b[j] = e[j] = b[j - 1] + c[j - 1];
 650         for (i = beg, j = 0; i != end;) {
 651                 rstype_t tmp = *i, swap;
 652                 int x;
 653                 for (;;) {
 654                         x = rskey(tmp)>>s&m;
 655                         if (e[x] == i) break;
 656                         swap = tmp; tmp = *e[x]; *e[x]++ = swap;
 657                 }
 658                 *i++ = tmp;
 659                 ++e[x];
 660                 while (j != size && i >= b[j]) ++j;
 661                 while (j != size && e[j-1] == b[j]) ++j;
 662                 if (i < e[j-1]) i = e[j-1];
 663         }
 664         if (s) {
 665                 s = s > n_bits? s - n_bits : 0;
 666                 for (j = 0; j < size; ++j) {
 667                         if (c[j] >= RS_MIN_SIZE) rs_sort2(b[j], e[j], n_bits, s);
 668                         else if (c[j] >= 2) rs_insertsort(b[j], e[j]);
 669                 }
 670         }
 671 }
 672 */
 673 void radix_sort(unsigned *array, int offset, int end, int shift) {
 674     int x, y, value, temp;
 675     int last[256] = { 0 }, pointer[256];
 676
 677     for (x=offset; x<end; ++x) {
 678         ++last[(array[x] >> shift) & 0xFF];
 679     }
 680
 681     last[0] += offset;
 682     pointer[0] = offset;
 683     for (x=1; x<256; ++x) {
 684         pointer[x] = last[x-1];
 685         last[x] += last[x-1];
 686     }
 687
 688     for (x=0; x<256; ++x) {
 689         while (pointer[x] != last[x]) {
 690             value = array[pointer[x]];
 691             y = (value >> shift) & 0xFF;
 692             while (x != y) {
 693                 temp = array[pointer[y]];
 694                 array[pointer[y]++] = value;
 695                 value = temp;
 696                 y = (value >> shift) & 0xFF;
 697             }
 698             array[pointer[x]++] = value;
 699         }
 700     }
 701
 702     if (shift > 0) {
 703         shift -= 8;
 704         for (x=0; x<256; ++x) {
 705             temp = x > 0 ? pointer[x] - pointer[x-1] : pointer[0] - offset;
 706             if (temp > 64) {
 707                 radix_sort(array, pointer[x] - temp, pointer[x], shift);
 708             } else if (temp > 1) rs_insertsort(array + pointer[x] - temp, array + pointer[x]);
 709         }
 710     }
 711 }
 712 /*************************
 713  *** END OF RADIX SORT ***
 714  *************************/
 715
 716 template< class _Type, unsigned long PowerOfTwoRadix, unsigned long Log2ofPowerOfTwoRadix, long Threshold >
 717 inline void _RadixSort_Unsigned_PowerOf2Radix_1( _Type* a, long last, _Type bitMask, unsigned long shiftRightAmount )
 718 {
 719         const unsigned long numberOfBins = PowerOfTwoRadix;
 720         unsigned long count[ numberOfBins ];
 721         for( unsigned long i = 0; i < numberOfBins; i++ )
 722                 count[ i ] = 0;
 723         for ( long _current = 0; _current <= last; _current++ ) // Scan the array and count the number of times each value appears
 724         {
 725                 unsigned long digit = (unsigned long)(( a[ _current ] & bitMask ) >> shiftRightAmount ); // extract the digit we are sorting based on
 726                 count[ digit ]++;
 727         }
 728         long startOfBin[ numberOfBins ], endOfBin[ numberOfBins ], nextBin;
 729         startOfBin[ 0 ] = endOfBin[ 0 ] = nextBin = 0;
 730         for( unsigned long i = 1; i < numberOfBins; i++ )
 731                 startOfBin[ i ] = endOfBin[ i ] = startOfBin[ i - 1 ] + count[ i - 1 ];
 732         for ( long _current = 0; _current <= last; )
 733         {
 734                 unsigned long digit;
 735                 _Type tmp = a[ _current ];  // get the compiler to recognize that a register can be used for the loop instead of a[_current] memory location
 736                 while ( true ) {
 737                         digit = (unsigned long)(( tmp & bitMask ) >> shiftRightAmount );   // extract the digit we are sorting based on
 738                         if ( endOfBin[ digit ] == _current )
 739                                 break;
 740                         _Type tmp2;
 741                         //_swap( tmp, a[ endOfBin[ digit ] ] );
 742                         tmp2 = a[endOfBin[digit]]; a[endOfBin[digit]] = tmp; tmp = tmp2;
 743                         endOfBin[ digit ]++;
 744                 }
 745                 a[ _current ] = tmp;
 746                 endOfBin[ digit ]++;   // leave the element at its location and grow the bin
 747                 _current++;  // advance the current pointer to the next element
 748                 while( _current >= startOfBin[ nextBin ] && nextBin < numberOfBins )
 749                         nextBin++;
 750                 while( endOfBin[ nextBin - 1 ] == startOfBin[ nextBin ] && nextBin < numberOfBins )
 751                         nextBin++;
 752                 if ( _current < endOfBin[ nextBin - 1 ] )
 753                         _current = endOfBin[ nextBin - 1 ];
 754         }
 755         bitMask >>= Log2ofPowerOfTwoRadix;
 756         if ( bitMask != 0 )   // end recursion when all the bits have been processes
 757         {
 758                 if ( shiftRightAmount >= Log2ofPowerOfTwoRadix ) shiftRightAmount -= Log2ofPowerOfTwoRadix;
 759                 else shiftRightAmount  = 0;
 760                 for( unsigned long i = 0; i < numberOfBins; i++ )
 761                 {
 762                         long numberOfElements = endOfBin[ i ] - startOfBin[ i ];
 763                         if ( numberOfElements >= Threshold )  // endOfBin actually points to one beyond the bin
 764                                 _RadixSort_Unsigned_PowerOf2Radix_1< _Type, PowerOfTwoRadix, Log2ofPowerOfTwoRadix, Threshold >( &a[ startOfBin[ i ]], numberOfElements - 1, bitMask, shiftRightAmount );
 765                         else if ( numberOfElements >= 2 )
 766                                 rs_insertsort(&a[ startOfBin[ i ]], &a[ endOfBin[ i ]]);
 767                 }
 768         }
 769 }
 770 inline void RadixSortInPlace_HybridUnsigned_Radix256( unsigned* a, unsigned long a_size )
 771 {
 772         if ( a_size < 2 ) return;
 773         unsigned long bitMask = 0xFF000000; // bitMask controls how many bits we process at a time
 774         unsigned long shiftRightAmount = 24;
 775         if ( a_size >= 32 )
 776                 _RadixSort_Unsigned_PowerOf2Radix_1<unsigned, 256, 8, 32>(a, a_size - 1, bitMask, shiftRightAmount );
 777         else
 778                 rs_insertsort(a, a + a_size);
 779 }
 780
 781 struct intcmp_t {
 782         inline int operator() (int a, int b) const {
 783                 return a < b? -1 : a > b? 1 : 0;
 784         }
 785 };
 786
 787 int compare_int(int a, int b)
 788 {
 789         return a < b? -1 : a > b? 1 : 0;
 790 }
 791 int compare(const void *a, const void *b)
 792 {
 793         return *((int*)a) - *((int*)b);
 794 }
 795
 796 int main(int argc, char *argv[])
 797 {
 798         int i, N = 50000000;
 799         int *array, *temp;
 800         clock_t t1, t2;
 801         if (argc == 1) fprintf(stderr, "Usage: %s [%d]\n", argv[0], N);
 802         if (argc > 1) N = atoi(argv[1]);
 803         temp = (int*)malloc(sizeof(int) * N);
 804         array = (int*)malloc(sizeof(int) * N);
 805
 806         srand48(11);
 807         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 808         t1 = clock();
 809         rs_sort((unsigned*)array, (unsigned*)array + N, 8, 24);
 810         t2 = clock();
 811         fprintf(stderr, "radix sort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 812         for (i = 0; i < N-1; ++i) {
 813                 if (array[i] > array[i+1]) {
 814                         fprintf(stderr, "Bug in radix sort!\n");
 815                         exit(1);
 816                 }
 817         }
 818         t1 = clock();
 819         rs_sort((unsigned*)array, (unsigned*)array + N, 8, 24);
 820         t2 = clock();
 821         fprintf(stderr, "radix sort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 822
 823         srand48(11);
 824         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 825         t1 = clock();
 826         RadixSortInPlace_HybridUnsigned_Radix256((unsigned*)array, N);
 827 //      radix_sort((unsigned*)array, 0, N, 24);
 828         t2 = clock();
 829         fprintf(stderr, "vd's radix sort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 830         for (i = 0; i < N-1; ++i) {
 831                 if (array[i] > array[i+1]) {
 832                         fprintf(stderr, "Bug in radix sort!\n");
 833                         exit(1);
 834                 }
 835         }
 836         t1 = clock();
 837         RadixSortInPlace_HybridUnsigned_Radix256((unsigned*)array, N);
 838 //      radix_sort((unsigned*)array, 0, N, 24);
 839         t2 = clock();
 840         fprintf(stderr, "vd's radix sort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 841
 842         srand48(11);
 843         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 844         t1 = clock();
 845         sort(array, array+N);
 846         t2 = clock();
 847         fprintf(stderr, "STL introsort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 848         t1 = clock();
 849         sort(array, array+N);
 850         t2 = clock();
 851         fprintf(stderr, "STL introsort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 852
 853         srand48(11);
 854         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 855         t1 = clock();
 856         stable_sort(array, array+N);
 857         t2 = clock();
 858         fprintf(stderr, "STL stablesort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 859         t1 = clock();
 860         stable_sort(array, array+N);
 861         t2 = clock();
 862         fprintf(stderr, "STL stablesort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 863
 864         srand48(11);
 865         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 866         t1 = clock();
 867         make_heap(array, array+N);
 868         sort_heap(array, array+N);
 869         t2 = clock();
 870         fprintf(stderr, "STL heapsort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 871         for (i = 0; i < N-1; ++i) {
 872                 if (array[i] > array[i+1]) {
 873                         fprintf(stderr, "Bug in heap_sort!\n");
 874                         exit(1);
 875                 }
 876         }
 877         t1 = clock();
 878         make_heap(array, array+N);
 879         sort_heap(array, array+N);
 880         t2 = clock();
 881         fprintf(stderr, "STL heapsort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 882
 883         srand48(11);
 884         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 885         t1 = clock();
 886         ks_combsort(int, N, array);
 887         t2 = clock();
 888         fprintf(stderr, "combsort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 889         for (i = 0; i < N-1; ++i) {
 890                 if (array[i] > array[i+1]) {
 891                         fprintf(stderr, "Bug in combsort!\n");
 892                         exit(1);
 893                 }
 894         }
 895
 896         srand48(11);
 897         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 898         t1 = clock();
 899         qsort(array, N, sizeof(int), compare);
 900         t2 = clock();
 901         fprintf(stderr, "libc qsort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 902
 903         srand48(11);
 904         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 905         t1 = clock();
 906         ks_introsort(int, N, array);
 907         t2 = clock();
 908         fprintf(stderr, "my introsort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 909         for (i = 0; i < N-1; ++i) {
 910                 if (array[i] > array[i+1]) {
 911                         fprintf(stderr, "Bug in intro_sort!\n");
 912                         exit(1);
 913                 }
 914         }
 915         t1 = clock();
 916         ks_introsort(int, N, array);
 917         t2 = clock();
 918         fprintf(stderr, "introsort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 919
 920         srand48(11);
 921         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 922         t1 = clock();
 923         ks_mergesort(int, N, array, 0);
 924         t2 = clock();
 925         fprintf(stderr, "iterative mergesort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 926         for (i = 0; i < N-1; ++i) {
 927                 if (array[i] > array[i+1]) {
 928                         fprintf(stderr, "Bug in merge_sort!\n");
 929                         exit(1);
 930                 }
 931         }
 932         t1 = clock();
 933         ks_mergesort(int, N, array, 0);
 934         t2 = clock();
 935         fprintf(stderr, "iterative mergesort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 936
 937         srand48(11);
 938         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 939         t1 = clock();
 940         ks_heapmake(int, N, array);
 941         ks_heapsort(int, N, array);
 942         t2 = clock();
 943         fprintf(stderr, "my heapsort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 944         for (i = 0; i < N-1; ++i) {
 945                 if (array[i] > array[i+1]) {
 946                         fprintf(stderr, "Bug in heap_sort!\n");
 947                         exit(1);
 948                 }
 949         }
 950         t1 = clock();
 951         ks_heapmake(int, N, array);
 952         ks_heapsort(int, N, array);
 953         t2 = clock();
 954         fprintf(stderr, "heapsort (sorted): %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 955
 956         srand48(11);
 957         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 958         t1 = clock();
 959         heapSort(array, N);
 960         t2 = clock();
 961         fprintf(stderr, "Paul's heapsort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 962         for (i = 0; i < N-1; ++i) {
 963                 if (array[i] > array[i+1]) {
 964                         fprintf(stderr, "Bug in intro_sort!\n");
 965                         exit(1);
 966                 }
 967         }
 968
 969         srand48(11);
 970         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 971         t1 = clock();
 972         quickSort(array, N);
 973         t2 = clock();
 974         fprintf(stderr, "Paul's quicksort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 975         for (i = 0; i < N-1; ++i) {
 976                 if (array[i] > array[i+1]) {
 977                         fprintf(stderr, "Bug in intro_sort!\n");
 978                         exit(1);
 979                 }
 980         }
 981
 982         srand48(11);
 983         for (i = 0; i < N; ++i) array[i] = (int)lrand48();
 984         t1 = clock();
 985         mergeSort(array, N);
 986         t2 = clock();
 987         fprintf(stderr, "Paul's mergesort: %.3lf\n", (double)(t2-t1)/CLOCKS_PER_SEC);
 988         for (i = 0; i < N-1; ++i) {
 989                 if (array[i] > array[i+1]) {
 990                         fprintf(stderr, "Bug in intro_sort!\n");
 991                         exit(1);
 992                 }
 993         }
 994
 995         free(array); free(temp);
 996         return 0;
 997 }