radix.hpp

#ifndef HAVE_RADIX_HPP
#define HAVE_RADIX_HPP
// Autogenerated - do not edit by hand !
#include <algorithm>
#include <vector>

namespace radix {

template <class T, class I>
struct radix {
  const std::vector<T>& x;
  std::vector<T> x_sort;
  std::vector<I> x_order;
  static const int radix_width = 8;
  static const int total_width = sizeof(T) * 8;
  static const int num_keys = (1 << radix_width);
  static const int mask = num_keys - 1;
  size_t key(T x, int k) { return (x >> k) & mask; }

  radix(const std::vector<T>& x) : x(x) {
    static_assert(T(-1) > T(0), "Value type must be unsigned");
    static_assert(total_width % radix_width == 0,
                  "Value width must be a multiple of radix width");
  }
  template <bool get_order>
  void run_sort() {
    T bitwise_min = ~0, bitwise_max = 0;
    for (size_t i = 0; i < x.size(); i++) {
      bitwise_min &= x[i];
      bitwise_max |= x[i];
    }

    x_sort = x;
    x_order.resize(x.size() * get_order);
    for (size_t i = 0; i < x_order.size(); i++) x_order[i] = i;

    std::vector<size_t> count(num_keys);
    std::vector<size_t> pos(num_keys);

    std::vector<I> y_order(x.size() * get_order);
    std::vector<T> y_sort(x.size());
    for (size_t k = 0; k < total_width; k += radix_width) {
      if (key(bitwise_min, k) == key(bitwise_max, k)) continue;

      std::fill(count.begin(), count.end(), 0);
      for (size_t i = 0; i < x.size(); i++) {
        count[key(x[i], k)]++;
      }

      std::fill(pos.begin(), pos.end(), 0);
      for (size_t i = 1; i < pos.size(); i++) {
        pos[i] = pos[i - 1] + count[i - 1];
      }
      for (size_t i = 0; i < x.size(); i++) {
        T x_sort_i = x_sort[i];
        size_t& j = pos[key(x_sort_i, k)];

        y_sort[j] = x_sort_i;
        if (get_order) y_order[j] = x_order[i];
        j++;
      }
      std::swap(x_sort, y_sort);
      std::swap(x_order, y_order);
    }
  }
  std::vector<T> sort() {
    run_sort<false>();
    return x_sort;
  }
  std::vector<I> order() {
    run_sort<true>();
    return x_order;
  }
  std::vector<I> first_occurance() {
    run_sort<true>();
    std::vector<I> ans(x_order.size());
    for (size_t i = 0; i < ans.size(); i++) ans[i] = i;
    for (size_t i = 1; i < x_sort.size(); i++) {
      if (x_sort[i - 1] == x_sort[i]) {
        ans[x_order[i]] = ans[x_order[i - 1]];
      }
    }
    return ans;
  }
};

template <class I, class T>
std::vector<I> order(const std::vector<T>& x) {
  return radix<T, I>(x).order();
}

template <class I, class T>
std::vector<I> first_occurance(const std::vector<T>& x) {
  return radix<T, I>(x).first_occurance();
}

template <class I, class T>
std::vector<I> factor(const std::vector<T>& x) {
  std::vector<I> y = first_occurance<I>(x);
  std::vector<I> tab(y.size());
  for (size_t i = 0, k = 0; i < y.size(); i++) {
    if (y[i] == i)
      tab[i] = k++;
    else
      tab[i] = tab[y[i]];
  }
  return tab;
}

}  // namespace radix
#endif  // HAVE_RADIX_HPP