#include <stdint.h>
#include <sys/time.h>
#include <cstdlib>
#include <stdexcept>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>

////////////////////////////////////////
//
// Returns current time in µsec
//
int64_t now() 
 {
  struct timeval tv;
  gettimeofday(&tv,0);
  return tv.tv_sec*((int64_t)1000000)+tv.tv_usec;
 }


template <typename T>
void cmpexchg(T & a, T & b) 
 { 
  using std::swap; // ADL safe
  if (a>b) swap(a,b); 
 }

template <typename T>
const T & median3( T & a, T & b, T & c )
 {
  cmpexchg(a,c);
  cmpexchg(a,b);
  cmpexchg(b,c);
  return b;
 }

template <typename T>
class med_heap
 {
 private:

  std::vector<T> & heap;

  int left_child(int i) { return 2*i+1; }
  int right_child(int i) { return 2*i+2; }
 
  void heapify()
   {
    for (int current=heap.size()/2-1;current>-1;current--)
     {
      int left=left_child(current);
      int right=right_child(current);

      // check if it has two children
      // (otherwise give up)
      //
      if ( (left < heap.size()) && 
           (right < heap.size()))
       {
        T & a = heap[current];
        T & b = heap[left];
        T & c = heap[right];

        const T & med = median3(a,b,c);

        if (b==med)
         std::swap(a,b);
        else
         if (c==med)
          std::swap(a,c);
         else
          ; // a is already the median
       }
      else
       ; // has only one child, so
         // already "median"
     }
   }

 public:

  size_t size() const { return heap.size(); }

  // lets you peek at the next
  const T & median() const
   {
    return heap[0];
   }


  med_heap(std::vector<T> & v)
   : heap(v)
   {
    heapify();
   }
 };

template <typename T>
void show( const med_heap<T> & this_heap)
 {
  int c=0;
  int med=this_heap.heap[0];
  for (int i=0;i<this_heap.heap.size();i++)
   {
    if (this_heap.heap[i]<=med) c++;
    std::cout << this_heap.heap[i] << " ";
   }
  std::cout
   << " ::: " << c
   << " ("
   << this_heap.size()/2
   << ")"
   << std::endl;

 }

////////////////////////////////////////
//
// Implements the generic selection algorithm
// loosely based on the quicksort partition
// algorithm (more or less any of the quicksort
// partition algorithms will do just fine)
//
// I'm sure there's a description of this
// algorithm or a variant of in Brassard and
// Bratley's introduction to algorithms (and
// it may be due to Wirth, from his book
// algorithms + data structures = programs)
//
// NOTE: array v gets modified in the process!
//
template <typename T>
const T & select( std::vector<T> & v,
                  int select_index)
  {
   // check if lo_index <= select_index <= hi_index
   // ...eventually
   //
   int lo_index=0;
   int hi_index=v.size()-1;

   while (hi_index>lo_index)
    {
     // classic mid-point
     // pivot selection (greatly
     // reduces chances of n^2
     // behavior if v is
     // nearly sorted, as it would
     // be after a few call to select.
     //
     const T pivot = v[(lo_index+hi_index)/2];

     int x=lo_index;
     int y=hi_index;

     // basic pivot/partition
     // algorithm
     //
     while (x<=y)
      {
       while (v[x]<pivot) x++;
       while (v[y]>pivot) y--;
       if (x<=y)
        {
         std::swap(v[x],v[y]);
         x++;
         y--;
        }
      }

     // check on which side of
     // the partition lands the
     // select_index, so that we
     // iterate only on the half
     // that contains the select_index
     //
     if (select_index >= x)
      lo_index=x;
     else
      hi_index=x-1;
    }

   // in our application, the actual
   // color isn't very important, unless
   // where at the last stage
   return v[select_index];
  }


////////////////////////////////////////
int main(int argc, char *argv[])
 {
  const size_t vector_size=1000;
  const int range=10000;


  for (int i=0;i<1000;i++)
   {
    std::vector<int> v(vector_size);

    for (int i=0;i<vector_size;i++)
     v[i]=std::rand() % range;

    std::vector<int> v2=v;


    int64_t start,stop;

    start=now();
    med_heap<int> z(v);
    int heap_med = z.median();
    stop=now();

    int64_t heap_med_time=stop-start;

    start=now();
    std::sort(v.begin(),v.end());
    int sort_med  = v[vector_size/2];
    stop=now();

    int64_t sort_med_time=stop-start;

    start=now();
    int select_med = select(v2,vector_size/2);
    stop=now();
    
    int64_t select_med_time=stop-start;

    std::cout
     << std::fixed
     << std::setprecision(7)
     << heap_med
     << " "
     << sort_med
     << " "
     << select_med
     << " "
     << (heap_med-sort_med)/(float)range
     << " "
     << heap_med_time/1000000.0
     << " "
     << sort_med_time/1000000.0
     << " "
     << select_med_time/1000000.0
     << std::endl;
   }

  return 0;
 }