Fast inplace sort of byte array

Question

I've got a small problem and can't find a satisfactory solution for it. There's a byte array and I need these bytes sorted by high 7 bits while preserving the order of low bits.

So originally it looked like this:

// sort buf[N] to tmp[N]
uint offs[128+1]; uint c,i,s;
for( i=0; i<128; i++ ) offs[i]=0;
for( i=0; i<l; i++ ) offs[buf[i]>>1]++;
for( i=0,s=0; i<128; i++ ) c=offs[i], offs[i]=s, s+=c; offs[i]=s;

byte* tmp = new byte[N];
for( i=0; i<N; i++ ) c=buf[i], tmp[offs[c>>1]++]=c; // sort

But these blocks are large enough (8M currently), and I want to use multiple threads, and an extra 8M per thread is noticeable.

So I tried to use some simple radix sort:

void radix( byte* buf, uint h, uint l, uint mask ) {
  uint p = (h+l)>>1, q = h; 
  uint i = offs[h], j = offs[l]-1; h = offs[p]; 
  if( (i<h) && (j>=h) ) {
    byte c = buf[i], d = buf[j];
    while( (i<h) && (j>=h) ) {
      while( (c&mask)==0 ) c = buf[++i]; // find value with bit 1
      while( (d&mask)!=0 ) d = buf[--j]; // find value with bit 0
      buf[i]=d; buf[j]=c; // swap 1-0 -> 0-1
      c = buf[++i]; d = buf[--j];
    }
    if( mask>=4 ) {
      radix( buf, q,p, mask>>1 );
      radix( buf, p,l, mask>>1 );
    }
  }
}

But it changes the order of these low bits and it becomes unusable.

Actually some simpler methods, like bubblesort, just do what I want, but they're much slower, and speed is an issue too.

So currently I sort smaller blocks via a temp buffer, then use an index table to access partially sorted chunks in order:

struct tmpsort {

  enum{ blocksize = (1<<16)-1 };

  unsigned short ofs[(max_quants+blocksize-1)/blocksize][probN];

  tmpsort( byte* buf, uint f_len ) {
    uint i,j,k;
    uint freq[2*probN]; // prob freqs
    byte tmp[blocksize+1];

    for( k=0,j=0; k<f_len; k+=blocksize,j++ ) {
      uint l = Min(k+blocksize,f_len)-k;
      byte* p = &buf[k];

      // compute offsets of sorted chunks
      for( i=0; i<2*probN; i++ ) freq[i]=0;
      for( i=0; i<l; i++ ) freq[p[i]]++;
      for( i=0; i<probN; i++ ) freq[i+1]=freq[2*i+0]+freq[2*i+1]; // 1=0+1, 2=2+3, 3=4+5
      freq[0] = 0;
      for( i=0; i<probN; i++ ) freq[i+1]+=freq[i];
      for( i=0; i<probN; i++ ) ofs[j][i]=freq[i+1];

      // sort the block via tmp
      for( i=0; i<l; i++ ) { byte c=p[i]; tmp[freq[c>>1]++]=c; }
      for( i=0; i<l; i++ ) p[i]=tmp[i];
    }
  }

};

[...]

tmpsort ts( buf, f_len );
for( i=0; i<probN; i++ ) {
  for( k=0,j=0; k<f_len; k+=ts.blocksize,j++ ) {
    uint x = i>0 ? ts.ofs[j][i-1] : 0;
    for(; x<ts.ofs[j][i]; x++ ) putc( buf[k+x],g );
  }
}

But tmp[] and ofs[] arrays use too much stack space, and its not a complete sort, so I keep wondering whether there's some neat solution for this.

A sample of data and my implementations are available here: http://nishi.dreamhosters.com/u/tmpsort_v0.rar

Answer 1

Why not just use any standard in-place, stable sorting algorithm, e.g. Insertion Sort, and implement an appropriate comparator function ?

Answer 2

This can be accomplished with relatively simple code in a little more than O(n log n) time using a version of radix sort that performs a stable sort on each of the 7 important bits, from least significant to most significant. The advantage of this technique relative to a stable in-place merge-sort is that the code is much simpler if you are writing it all yourself.

Here is the function to perform an in-place stable sort by one specified bit. Here, it is written recursively for simplicity using O(lg n) stack space (this stack space usage can be eliminated if you want by using a for loop to organize the divide and conquer approach):

// sort array x from i to j by bit b
sort(x, i, j, b) {
  if (i >= j - 1) return;
  mid = (i + j) / 2;
  sort(x, i, mid, b);
  sort(x, mid, j, b);
  first1 = -1;
  last0 = -1;
  for (k = i; k < j; k++) {
    if (first1 < 0 && isSet(x[k], b)) first1 = k;
    if (!isSet(x[k], b)) last0 = k;
  }
  if (last0 < first1) return;

  // the sequence of bit b generally looks something like 0000011100000111111
  // so we reverse from the first 1 to the last 0
  reverse(x, first1, last0afterfirst1);
  newlast0 = first1;
  while (!isSet(x[++newlast0], b));
  newlast0--;

  // the elements in the range first1..last0 are in the wrong order, so reverse
  reverse(x, first1, newlast0);
  reverse(x, newlast0 + 1, last0);
}

The function isSet tests whether a bit is set and reverse performs in-place array reversal. The above sorting subroutine is called on each bit as follows (as in radix sort):

sort(x) {
  for (b = 1; b < 8; b++) {
    sort(x, 0, n, b);
  }
}

The total running time is "O(7 * n log n)". The extra factor of 7 could be variable if this algorithm were generalized.