radixsort.cpp

template <class IntegralType>
inline size_t idx(IntegralType x, size_t p, size_t d) {
  // e.g. p=1 d=2
  // 101101 -> 11
  return ((x >> (p * d))) & ((1 << d) - 1);
}

template <class BidirectionalIterator>
void radix_sort(BidirectionalIterator first, BidirectionalIterator last,
                int d = 8, int w = 64) {
  using ValueType =
      typename std::iterator_traits<BidirectionalIterator>::value_type;
  static_assert(std::is_integral<ValueType>::value);

  std::vector<size_t> count(1 << d);
  std::vector<ValueType> a(first, last);
  std::vector<ValueType> tmp(a.size());

  for (size_t p = 0; p < w / d; p++) {
    std::fill(count.begin(), count.end(), 0);
    for (auto x : a) ++count[idx(x, p, d)];
    std::partial_sum(count.begin(), count.end(), count.begin());
    for (auto it = a.rbegin(); it != a.rend(); ++it) {
      tmp[--count[idx(*it, p, d)]] = *it;
    }
    std::swap(a, tmp);
  }
  std::copy(a.begin(), a.end(), first);
}