Fast inplace sort of byte array

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<开发者_如何学Gol; 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

---

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

---

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.

---

Having additional 64kB, you can (as you have noticed) store a 512 kbit block (minus some fixed amount of indexing data) in compressed form (storing only the lowest bits for every key) Go over the large blocks and convert them to their compressed-sorted forms, compacting them as you go at the beginning of the whole array.

Now merge the compressed forms into one big compressed form (easy with the 7M freed up.) Then uncompress back to the sorted array.

This is O(N), although the constant looks rather large with 3 passes that involve some nontrivial bit operations.

---

It's possible to implement quicksort as a stable sort. In terms of big-O, it's no better than insertion sort, but in practice it will perform a lot better. If you hard-code sorting networks for the leaves of size up to 6 or 8, I think that's about the best performance you're going to get for a stable, in-place sort.

Actually... supposedly there's such a thing as an in-place, stable merge sort. In terms of ideal theoretical characteristics, it's the holy grail of sorting - in-place, true O(n log n), and stable, all at the same time. But I suspect it's a huge pain to implement and has rather large constant terms to go with that big-O.