Q:

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<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

A:

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.

