trefny.net

rsort.h

#if !defined(CGL_RSORT_H)
# define CGL_RSORT_H

#include "config.h"

#include <vector>

namespace cgl {

class rsort
{
 public :

  rsort & sort(const int   * data_p, unsigned int size);
  rsort & sort(const float * data_p, unsigned int size);

  const std::vector<unsigned int> & pos() const { return positions; }

 private :

  unsigned int counters[4 * 256];
  unsigned int offsets [256];

  std::vector<unsigned int> positions;
};

} // cgl

rsort.cpp

#include "rsort.h"

#include <assert.h>

cgl::rsort &
cgl::rsort::sort(const int * data_p, unsigned int size)
{
  assert(data_p && size); // FAILURE: null input

  // clear counters
  for(unsigned int i = 0; i < 4 * 256; ++i) counters[i] = 0;

  // init positions
  if(positions.size() < size)
    for(unsigned int i = 0; i < size; ++i) positions.push_back(i);

  bool already_sorted = true;
  int  negvals        = 0;
  for(unsigned int i = 0; i < size; ++i)
  {

    // check for sorted input
    if(i && data_p[i] < data_p[i-1])
      already_sorted = false;

    // check for presence of negative values
    if(data_p[i] < 0)
      negvals++;

    // init counters for all 4 passes
    ++counters[0   + ((data_p[i] & 0x000000ff) >> 0)];
    ++counters[256 + ((data_p[i] & 0x0000ff00) >> 8)];
    ++counters[512 + ((data_p[i] & 0x00ff0000) >> 16)];
    ++counters[768 + ((data_p[i] & 0xff000000) >> 24)];
  }

  if(already_sorted)
    return *this;

  std::vector<unsigned int> new_positions(positions.size());
  for(int pass = 0; pass < 4; ++pass)
  {
    if(!negvals || pass < 3)
    {
      // compute offsets
      offsets[0] = 0;
      for(int i = 1; i < 256; ++i)
        offsets[i] = offsets[i - 1] + counters[(pass << 8) + i - 1];
    }
    else
    {
      // positive values' offsets are biased to be sorted
      // in positions after negative values
      offsets[0] = negvals;
      for(int i = 1; i < 128; ++i)
        offsets[i] = offsets[i - 1] + counters[(pass << 8) + i - 1];
      offsets[128] = 0;
      for(int i = 129; i < 256; ++i)
        offsets[i] = offsets[i - 1] + counters[(pass << 8) + i - 1];
    }

    // sort positions
    for(unsigned int i = 0; i < size; ++i)
    {
      unsigned int pos  = positions[i];
      unsigned int bval = (data_p[pos] >> (8 * pass)) & 0xff;
      new_positions[offsets[bval]++] = pos;
    }

    // save new positions for the next pass
    positions.swap(new_positions);
  }

  return *this;
}

// IEEE floats are required for the following to work
// an IEEE float is a 32-bit number with the following structure:
//        [sign-1bit|exponent-8bits|mantissa-23bits]
//  sign bit is set for negative values; exponent is biased and
//  thus can be treated as always positive
cgl::rsort &
cgl::rsort::sort(const float * data_p, unsigned int size)
{
  assert(data_p && size); // FAILURE: null input

  // clear counters
  for(unsigned int i = 0; i < 4 * 256; ++i) counters[i] = 0;

  // init positions
  if(positions.size() < size)
    for(unsigned int i = 0; i < size; ++i) positions.push_back(i);

  bool already_sorted = true;
  int  negvals        = 0;
  for(unsigned int i = 0; i < size; ++i)
  {
    // check for sorted input
    if(i && data_p[i] < data_p[i-1])
      already_sorted = false;

    // check for presence of negative values
    if(data_p[i] < 0)
      negvals++;

    // init counters for all 4 passes
    unsigned int f2ui = *reinterpret_cast(&data_p[i]);
    ++counters[0   + ((f2ui & 0x000000ff) >> 0)];
    ++counters[256 + ((f2ui & 0x0000ff00) >> 8)];
    ++counters[512 + ((f2ui & 0x00ff0000) >> 16)];
    ++counters[768 + ((f2ui & 0xff000000) >> 24)];
  }

  if(already_sorted)
    return *this;

  std::vector<unsigned int> new_positions(positions.size());
  for(int pass = 0; pass < 4; ++pass)
  {
    if(pass < 3)
    {
      // compute offsets
      offsets[0] = 0;
      for(int i = 1; i < 256; ++i)
        offsets[i] = offsets[i - 1] + counters[(pass << 8) + i - 1];

      // sort positions
      for(unsigned int i = 0; i < size; ++i)
      {
        unsigned int pos  = positions[i];
        unsigned int bval = (*reinterpret_cast(&data_p[pos]) >> (8 * pass)) & 0xff;
        new_positions[offsets[bval]++] = pos;
      }
    }
    else
    {
      // positive values' offsets are biased to be
      // sorted in positions after negative values
      offsets[0] = negvals;
      for(int i = 1; i < 128; ++i)
        offsets[i] = offsets[i - 1] + counters[768 + i - 1];
      // negative values must have a reversed sorting order
      // (higher the exponent lower the position)
      offsets[255] = 0;
      for(int i = 254; i >= 128; --i)
        offsets[i] = offsets[i + 1] + counters[768 + i + 1];
      // radix sort is stable
      for(int i = 128; i < 256; ++i)
        offsets[i] += counters[768 + i];

      // sort positions
      for(unsigned int i = 0; i < size; ++i)
      {
        unsigned int pos  = positions[i];
        unsigned int bval = *reinterpret_cast(&data_p[pos]) >> 24;

        if(bval < 128) // positive/negative
          new_positions[offsets[bval]++] = pos;
        else
          new_positions[--offsets[bval]] = pos;
      }
    }

    // save new positions for the next pass
    positions.swap(new_positions);
  }

  return *this;
}