/*
  This is a one-pass version of the selection algorithm (LazySelect)
  described in [Rajeev Motwani and Prabhakar Raghavan, Randomized
  Algorithms, Cambridge University Press (1995)].

  Authors: Jesper Bojesen, Jyrki Katajainen, Lars Yde © 1999, 2002, 2011
*/

#ifndef __CPHSTL_RANDOMIZED_POSITIONER__
#define __CPHSTL_RANDOMIZED_POSITIONER__

#include <algorithm> // std::min, std::max
#include <cmath> // ilogb
#include <cstdlib> // std::rand
#include <functional> // std::binary_function
#include <iterator> // std::iterator_traits 
#include <vector> // std::vector
#include <utility>  // std::pair

extern int ilogb(double) throw();

namespace cphstl {

  namespace {
    
    template <typename R, typename E, typename integer>
    void tagged_sampler(R first, R past, 
			std::vector<std::pair<E, R> >& sample, integer n) {
      assert(sample.empty());
      typename std::vector<std::pair<E, R> >::iterator i;
      integer N = past - first;
      assert(RAND_MAX >= N);
      for (integer i = 0; i != n; ++i) {
	integer rnd = std::rand() % N;
	sample.push_back(std::pair<E, R>(*(first + rnd), first + rnd));
      }
    }

    template <typename E, typename R, typename C>
    class tagged_relation 
      : public std::binary_function<std::pair<E, R>, std::pair<E, R>, bool> {
    public:
      
      typedef std::pair<E, R> argument_type;
      
      explicit tagged_relation(C const& f) : less(f) {
      }
      
      bool operator() (argument_type const& x, argument_type const& y) const { 
	return less(x.first, y.first)/* ||
                                        (!less(y.first, x.first) && x.second < y.second)*/;
      }
      
    protected:
      
      C less;
    };

    template <typename E, typename R, typename I, typename C>
    std::pair<E, R> maximum(I first, I past, C less) {
      assert(first != past);
      E max = (*first).first;
      I location = first;
      I q = first;
      for (++q; q != past; ++q) {
	E x = (*q).first;
        if (! less(x, max)) {
	  max = x;
	  location = q;
        }
      }
      return *location;
    }

    template <typename E, typename R, typename I, typename C>
    std::pair<E, R> minimum(I first, I past, C less) {
      assert(first != past);
      E min = (*first).first;
      I location = first;
      I q = first;
      for (++q; q != past; ++q) {
	E x = (*q).first;
        if (less(x, min)) {
	  min = x;
	  location = q;
        }
      }
      return *location;
    }

    template <typename R, typename C>
    void two_way_partitioner(R first, R past, R p, C const& less) {
      typename std::iterator_traits<R>::value_type piv(*p);
      std::swap(*p, *(past-1));
      R store = first;
      for(R it=first; it!=past-1; ++it) {
          if (less(*it, piv)) {
              std::swap(*store, *it);
              ++store;
          }
      }
      std::swap(*(past-1), *store);
    }

    template <typename R>
    void rotate(R a, R b, R c) { 
      typedef typename std::iterator_traits<R>::value_type E;

      E tmp = *a;
      *a = *b; 
      *b = *c; 
      *c = tmp;
    }
  
    template <typename E, typename R, typename C>
    std::pair<R, R> three_way_partitioner(R first, R past, 
      std::pair<E, R> pivot1, std::pair<E, R> pivot2, C const& less) {
      R less_p1 = first;
      R less_p2 = first;
      R greater_p2 = past - 1;
      R greater_p1 = past - 1;
      
      while (true) {
	while (less(std::pair<E, R>(*less_p1, less_p1), pivot1)) {
	  less_p1++;
	}
	while (less(pivot2, std::pair<E, R>(*greater_p2, greater_p2))) {
	  greater_p2--;
	}
	if (less_p1 >= greater_p2) {
	  break;
	}
	if (less(pivot2, std::pair<E, R>(*less_p1, less_p1))) {
	  if (less(std::pair<E, R>(*greater_p2, greater_p2), pivot1)) {
	    std::iter_swap(less_p1++, greater_p2--);
	  }
	  else {
	    rotate(greater_p2--, less_p1++, less_p2++);
	  }
	}
	else {
	  if (less(std::pair<E, R>(*greater_p2, greater_p2), pivot1)) {
	    rotate(less_p1++, greater_p2--, greater_p1--);
	  }
	  else {
	    std::iter_swap(greater_p1--, greater_p2--);
	    std::iter_swap(less_p1++, less_p2++);
	  }
	}
      }
      if (less_p1 == greater_p2) {
	less_p1++; 
	greater_p2--;
      };
      std::swap_ranges(first, 
		       std::min(less_p2, greater_p2 + 1 - (less_p2 - first)), 
		       std::max(less_p2, greater_p2 + 1 - (less_p2 - first)));
      std::swap_ranges(less_p1,
		std::min(greater_p1 + 1, past - (greater_p1 + 1 - less_p1)),
		std::max(greater_p1 + 1, past - (greater_p1 + 1 - less_p1)));
      return std::make_pair(greater_p2 + 1 - (less_p2 - first),
			    past - (greater_p1 + 1 - less_p1));
    } 

    template <typename I>
    void show(I first, I past) {
      for (I p = first; p != past; ++p) {
	std::cout << *p << " ";
      }    
      std::cout << std::endl;
    }

  }

  template <typename R, typename integer, typename C>
  void randomized_positioner(R first, R past, integer k, C const& less) {
    typedef typename std::iterator_traits<R>::value_type E;
    typedef std::pair<E, R> tagged_element;

    tagged_relation<E, R, C> tagged_less = tagged_relation<E, R, C>(less);
    integer n = past - first;
    if (n < 8) {
      std::sort(first, past, less);
      return;
    }
    integer fourth_root = powl(n,1.0/4.0);
    integer square_root = powl(n,1.0/2.0);
    integer sample_size = powl(n,3.0/4.0);
    std::vector<tagged_element> sample;
    tagged_sampler(first, past, sample, sample_size);

    if (k <= fourth_root) {
      tagged_element p = maximum<E, R>(sample.begin(), sample.end(), less);
      R middle = first + sample_size;
      --middle;
      std::swap(*p.second, *middle);
      two_way_partitioner(first, past, middle, less);
      std::sort(first, middle, less); 
    }
    else if (k >= (n - fourth_root)) {
      tagged_element p = minimum<E, R>(sample.begin(), sample.end(), less);
      R middle = past - sample_size;
      std::swap(*p.second, *middle);
      two_way_partitioner(first, past, middle, less);
      std::sort(middle, past, less); 
    }
    else {
      integer x = k / fourth_root;
      integer ell = std::max(1, x - square_root);
      integer hoo = std::min(x + square_root, sample_size-1);
      typedef typename std::vector<tagged_element>::iterator I;
      std::sort(sample.begin(), sample.end(), tagged_less);
      tagged_element a = sample[ell];
      tagged_element b = sample[hoo];
      std::pair<R, R> m = three_way_partitioner(first, past, a, b, tagged_less);

      if (k < (m.first - first)) {
	std::sort(first, m.first, less);
      }
      else if (k > (m.second - first)) {
	std::sort(m.second, past, less);
      }
      else {
	std::sort(m.first, m.second, less);
      }
    }
  }
}

#endif