Radix sort:
===========
recursive implementation:

RRS( a, b, e, pos ) {
    if ( pos == -1 ) return; // all bits used

    q = partition_by_bit( a, b, e, pos ); // moves all elements with ZERO bit at position pos to the left of elements with same bit set to ONE
                                          // returns index of the first element with bit ONE

    RRS( a, b, q, pos-1 );
    RRS( a, q, e, pos-1 );
}


RRS( a, 0, size, sizeof(int)*8 -1 )

101 111 011 001 110 100
to start RRS( a, 0, 6, 3-1 )


011 001 - 101 111 110 100 // split by 2nd bit
-LEFT-    -----RIGHT----

001 - 011 //split LEFT by 1st bit 
-L1-  -L2-
and
100 101 - 111 110  //split RIGHT by 1st bit 
--R1--    --R2--    

L1: 001 - q=e 
L2: 011 - q=b
R1: 100 - 101 //split by 0th bit
R2: 110 - 111 //split by 0th bit


Notice that the structure of the algorithm is identical to the quicksort. We know how to parallelize it.


Straight radix sort:
====================

using digits as blocks:
1) looking at the last digit
16 82 89 63 16a 79 72 75 19 44
                                        19
            72              16a         79
            82  63  44  75  16          89
bins 0   1   2   3   4   5   6   7   8   9
     0   0   2   1   1   1   2   0   0   3
	 0   0   2   3   4   5   7   7   7   10
 
read from bins in bottom up order
82 72 63 44 75 16 16a 89 79 19

2) looking at the next digit (from last)
         19                      79
         16a                     75  89
         16  26      44      63  72  82
bins 0   1   2   3   4   5   6   7   8   9

read from bins in bottom up order
16 16a 19 26 44 63 72 75 79 82 89

DONE.

The order is stable - if two elements are equal,
     they will be in the same order in sorted array, as in the original.

Typically one will use blocks of bits instead of digits (faster):

////////////////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <algorithm>
#include <iostream>
#include <cassert>

int bits( int v, int pos, int len ) {
    return (v >> pos) % (1<<len);
}

int const width = 4;
int const num_bins     = 1<<width;

void StraightRadixSort(int *a, int size, int b)
{
    int count[num_bins];
    int *ta = new int[size];

    std::cerr << "b / width = " << b / width << std::endl;
    
    for ( int pass = 0; pass < (b / width); ++pass ) {
        for (int i=0;i<size;++i) { std::cout << a[i] << " "; } std::cerr << std::endl;

        for ( int j = 0; j < num_bins; ++j ) count[j] = 0;
       
        for ( int i = 0; i < size; ++i ) {
            std::cerr << "bits("<<a[i]<<","<<pass<<"*"<<width<<", "<<width<<" ) = " << bits(a[i],pass*width, width ) << std::endl;
            int bin_index = bits(a[i], pass*width, width);
            ++count[ bin_index ];
        }
        --count[0]; // 0-based indexing in target array 

        std::cerr << "            : "; for ( int j = 0; j < num_bins; ++j ) { std::cerr << std::setw(3) << j << " "; } std::cerr << std::endl;
        std::cerr << "counts      : "; for ( int j = 0; j < num_bins; ++j ) { std::cerr << std::setw(3) << count[j] << " "; } std::cerr << std::endl;
        for ( int j = 1; j < num_bins; j++) { count[j] += count[j-1]; }  // prefix sum on counts
        std::cerr << "prefix scan : "; for ( int j = 0; j < num_bins; ++j ) { std::cerr << std::setw(3) << count[j] << " "; } std::cerr << std::endl;
        
        // counts - after *inclusive* prefix scan gives the *last* position of the group
        // therefore going backwards
        for ( int i = size-1; i >= 0; --i) {
            int bin_index = bits(a[i], pass*width, width);
//            std::cerr << "ta[ "<<count[bin_index]<<" ] = "<<a[i] << std::endl;
            ta[ count[bin_index] ] = a[i];
            --count[bin_index];
        }
        
        // swap pointers to swap arrays to sort in next pass
        std::swap( ta, a );
    }
    //free(ta); // only correct if number of passes is even, needs a copy back to argument "a" if odd 
}

int main () {
    int a[8] = {805,306,958,900,937,715,304,992};
//    for (int i=0;i<100;++i) { a[i]=i; }
//    std::random_shuffle( a, a+100 );
    StraightRadixSort( a, sizeof(a)/sizeof(*a), sizeof(int)*8 );
    for (int i=0;i<sizeof(a)/sizeof(*a);++i) { std::cout << a[i] << " "; }
    std::cout << std::endl;
}
////////////////////////////////////////////////////////////////////////////////


805 306 958 900 937 715 304 992 
bits(805,0*4, 4 ) = 5
bits(306,0*4, 4 ) = 2
bits(958,0*4, 4 ) = 14
bits(900,0*4, 4 ) = 4
bits(937,0*4, 4 ) = 9
bits(715,0*4, 4 ) = 11
bits(304,0*4, 4 ) = 0
bits(992,0*4, 4 ) = 0
counts      : 2 0 1 0 1 1 0 0 0 1 0 1 0 0 1 0 
after --counts[0]
counts      : 1 0 1 0 1 1 0 0 0 1 0 1 0 0 1 0 

bit pattern   0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
prefix scan : 1 1 2 2 3 4 4 4 4 5  5  6  6  6  7  7 

Example: 992, bit pattern 0, so position 1 = counts[1]
 0  1  2  3  4  5  6  7  
 - 902 -  -  -  -  -  -  -  
 decrement corresponding count  --counts[1]

bit pattern   0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
prefix scan : 0 1 2 2 3 4 4 4 4 5  5  6  6  6  7  7 

next 304, bit pattern 0, so position 0 (was just updated!)
 0   1  2  3  4  5  6  7  
304 902 -  -  -  -  -  - 

Note that scanning backwards preserves stability.

It is possible to modify the way prefix scan is applied
so that arrays is scanned from left to right.

.....

After first pass:
304 992 306 900 805 937 715 958 

bits(304,1*4, 4 ) = 3
bits(992,1*4, 4 ) = 14
bits(306,1*4, 4 ) = 3
bits(900,1*4, 4 ) = 8
bits(805,1*4, 4 ) = 2
bits(937,1*4, 4 ) = 10
bits(715,1*4, 4 ) = 12
bits(958,1*4, 4 ) = 11
counts      : -1 0 1 2 0 0 0 0 1 0 1 1 1 0 1 0 
prefix scan : -1 -1 0 2 2 2 2 2 3 3 4 5 6 6 7 7 

position with -1 will never be used (they are -1 since there are no array elements 
with such bit pattern)

805 304 306 900 937 958 715 992 
bits(805,2*4, 4 ) = 3
bits(304,2*4, 4 ) = 1
bits(306,2*4, 4 ) = 1
bits(900,2*4, 4 ) = 3
bits(937,2*4, 4 ) = 3
bits(958,2*4, 4 ) = 3
bits(715,2*4, 4 ) = 2
bits(992,2*4, 4 ) = 3
counts      : -1 2 1 5 0 0 0 0 0 0 0 0 0 0 0 0 
prefix scan : -1 1 2 7 7 7 7 7 7 7 7 7 7 7 7 7 
304 306 715 805 900 937 958 992 
bits(304,3*4, 4 ) = 0
bits(306,3*4, 4 ) = 0
bits(715,3*4, 4 ) = 0
bits(805,3*4, 4 ) = 0
bits(900,3*4, 4 ) = 0
bits(937,3*4, 4 ) = 0
bits(958,3*4, 4 ) = 0
bits(992,3*4, 4 ) = 0
counts      : 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
prefix scan : 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
304 306 715 805 900 937 958 992 
bits(304,4*4, 4 ) = 0
bits(306,4*4, 4 ) = 0
bits(715,4*4, 4 ) = 0
bits(805,4*4, 4 ) = 0
bits(900,4*4, 4 ) = 0
bits(937,4*4, 4 ) = 0
bits(958,4*4, 4 ) = 0
bits(992,4*4, 4 ) = 0
counts      : 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
prefix scan : 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
304 306 715 805 900 937 958 992 
bits(304,5*4, 4 ) = 0
bits(306,5*4, 4 ) = 0
bits(715,5*4, 4 ) = 0
bits(805,5*4, 4 ) = 0
bits(900,5*4, 4 ) = 0
bits(937,5*4, 4 ) = 0
bits(958,5*4, 4 ) = 0
bits(992,5*4, 4 ) = 0
counts      : 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
prefix scan : 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
304 306 715 805 900 937 958 992 
bits(304,6*4, 4 ) = 0
bits(306,6*4, 4 ) = 0
bits(715,6*4, 4 ) = 0
bits(805,6*4, 4 ) = 0
bits(900,6*4, 4 ) = 0
bits(937,6*4, 4 ) = 0
bits(958,6*4, 4 ) = 0
bits(992,6*4, 4 ) = 0
counts      : 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
prefix scan : 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
304 306 715 805 900 937 958 992 
bits(304,7*4, 4 ) = 0
bits(306,7*4, 4 ) = 0
bits(715,7*4, 4 ) = 0
bits(805,7*4, 4 ) = 0
bits(900,7*4, 4 ) = 0
bits(937,7*4, 4 ) = 0
bits(958,7*4, 4 ) = 0
bits(992,7*4, 4 ) = 0
counts      : 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
prefix scan : 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
304 306 715 805 900 937 958 992 

To create a parallel version of this algorithm:

for ( passes ) {
    calculate bins // array may be split in any number of chunks and threads count in parallel, assuming thread-safe container for counts
    prefix scan    // parallel version covered earlier in this course
    scan original array and use values in bins for positions // array may be split in any number of chunks and threads copy elements in parallel, assuming thread-safe container for counts
}

there is a problem though:
we may (most probably will) loose stability of the sort
if thread T4 finds element with bit-pattern (bp) earlier then T1 (becasue it is closer to the starting end of the segment/chunk), 
then T4 will place found element into final array after the one found by T1 (remember scan is going backwards).


assume chunk for T1 is to the left of T4's chunk in the original array 
assume E1 and E2 are equal, in stable sort E1 will always stay to the left of E2, but:
time 0
chunk for T1                    chunk for T4
---------E1--                   ----E2----------
           <<                                 <<  direction of scan
           T1                                 T4
time 1
chunk for T1                    chunk for T4
---------E1--                   ----E2----------
         <<                                 <<  direction of scan
         T1                                 T4
T1 found E, and copies it to final array

final array
------------------------------------E1-------
               { block for element E }
time 2
chunk for T1                    chunk for T4
---------E1--                   ----E2----------
 <<                                 <<  direction of scan
         T1                                 T4
                                            T4 found E, and copies it to final array

final array
----------------------------------E2E1-------
               { block for element E }
NOT STABLE.


To solve the problem.
Solution 1) Keep the third line (in pseudo code) single threaded -- obviously correct, but slow.

Solution 2) Keep track of ranks ( Not really a solution, but there is an idea that may be developed into a solution:)
original array
     A1 C A2 B A3
T1   1  0 1  0 1   (looks for A's)
T2   0  0 0  1 0   (looks for B's)
T3   0  1 0  0 0   (looks for C's)
prefix scan 
     A1 C A2 B A3
T1   1  0 2  0 3   (looks for A's) <-- gives us the info how to order A's (A3 is the third of all A's)
T2   0  0 0  1 0   (looks for B's)
T3   0  1 0  0 0   (looks for C's)

so that if for final copy implemenetation we use 2 threads, and the first thread finds A2 before secondsees A3,
the first thread calculates proper position for A2:
    position (from the single threaded) -(total count of As) + rank(A2) = pos -3 + 2 = pos -1
    DO NOT decrement position!
then second thread finds A3 and
    position (from the single threaded) -(total count of As) + rank(A3) = pos -3 + 3 = pos 

so that 
   pos-1  pos
   A2     A3

The code works but is very demanding:
1) need an array of the same size as the original size for each character - absolutely unreasonable
2) each thread traverses the whole array ( may be parallelized thought)
3) number of threads is limited by the size of the alphabet



Another solution (better):
==========================
1) threads work on chunks
2) number of threads is any
3) additional space is the <number of bins> * <number of threads>

Here is an example:

13 09 95 84 71 | 29 64 80 05 60 | 91 29 76 37 97 | 26 52 87 14 84 <- array
 3  9  5  4  1 |  9  4  0  5  0 |  1  9  6  7  7 |  6  2  7  4  4 <- last digit
 

each thread uses its own array of bins:

per thread bins
T1:   0   1   0   1   1   1   0   0   0   1 
T2:   2   0   0   0   1   1   0   0   0   1 
T3:   1   0   0   0   0   0   1   2   0   1 
T4:   0   0   1   0   2   0   1   1   0   0 

Gcounts ( Gcounts[i] = T1[i]+T2[i]+T3[i]+T4[i] )
index 0   1   2   3   4   5   6   7   8   9
      3   1   1   1   4   2   2   3   0   3 

Fcounts ( exclusive scan of Gcounts)
index 0   1   2   3   4   5   6   7   8   9
      0   3   4   5   6  10  12  14  17  17 
Fcounts[i] is the first index for the group of element that match bit pattern "i":

Now each thread adjusts its bins from the first step by looking
at Gcounts and bins from the previous threads:
    for ( int i=0; i<10; ++i ) { // single threaded solution needs "if"s
        if ( T4[i] != 0 ) { T4[i] = fc[i]+T1[i]+T2[i]+T3[i]; }
        if ( T3[i] != 0 ) { T3[i] = fc[i]+T1[i]+T2[i]; }
        if ( T2[i] != 0 ) { T2[i] = fc[i]+T1[i]; }
        if ( T1[i] != 0 ) { T1[i] = fc[i]; }
    }
index 0   1   2   3   4   5   6   7   8   9
T1:   0   3   0   5   6  10   0   0   0  17 
T2:   0   0   0   0   7  11   0   0   0  18 
T3:   2   0   0   0   0   0  12  14   0  19 
T4:   0   0   4   0   8   0  13  16   0   0 
if each thread is doing its own part then the "if"s are not needed - those 
positions will never be looked at

examples: 
T1 updates position 4 as follows:
T1[4] = Fcounts[4]  = 6
so that element ending with 4 found by T1 will be placed at position 6 (there will 1 such element)

T2 updates position 4 as follows:
T2[4] = Fcounts[4]  + T1[4] = 6+1 = 7
so that element ending with 4 found by T2 will be placed at position 7 (there will 1 such element)

T3 skips position 4:

T4 updates position 4 as follows:
T4[4] = Fcounts[4] + T1[4] + T2[4] + T3[4] + T4[4] = 6+1+1+0 = 8
so that element ending with 4 found by T4 will be placed at position 8 (there will 2 such elements)

Note that stability requirement is satisfied, since threads are given chunks from left to right,
and T2 will be placing 64 after T1's 84. Also 14 and 84 found by T4 will be in the right order
since T4 scan its chunk left to right.