#ifndef D2SORTVEC_H
#define D2SORTVEC_H

#include <cassert>
#include <algorithm>
#include <cmath>
using namespace std;

#include <functionalobjectoperators.h>
#include <typedefs.h>

/*
\brief Sort a 2D vector.

*/
template< typename T, typename CX, typename CY >
class d2sortvec
{
  CX cmpx;
  CY cmpy;

  /** width:=N^(1/2). */
  void setwidth2D();
    
  /** width:=N^(2/3). */
  void setwidth3D();
public:

  /** The vector. */
  T * vi;
  /** The size of the vector. */
  uint N;

  /** The number of points within an interval. */
  uint width;
  /** Successive interval have opposite parity.
      So the snake top down, then botto up. */
  uint parity;

  d2sortvec(T * _vi, uintc _N, uintc _parity=0);

  d2sortvec
  (
    T * _vi, 
    uintc _N, 
    CX const & _cmpx, 
    CY const & _cmpy, 
    uintc _parity=0
  );

  /** Two dimensional order of n elements. */
  void eval();

  /** Three dimensional order using two dimensional 
      blocks. */
  void eval(uintc n2);

};

template< typename T, typename CX, typename CY >
void d2sortvecfunction
(
  T* vi,
  uintc N,
  CX const & cx,
  CY const & cy
)
{
  
  d2sortvec<T,CX,CY> s(vi,N,cx,cy);
  //s.eval();
}
    
  


//----------------------------------------------------------
//  Implementation



template< typename T, typename CX, typename CY >
d2sortvec<T,CX,CY>::d2sortvec
(
  T * _vi, 
  uintc _N, 
  uintc _parity
)
  : vi(_vi), N(_N), parity(_parity) 
{
  setwidth2D();
}

template< typename T, typename CX, typename CY >
d2sortvec<T,CX,CY>::d2sortvec
(
  T * _vi, 
  uintc _N, 
  CX const & _cmpx, 
  CY const & _cmpy, 
  uintc _parity
)
  : cmpx(_cmpx), cmpy(_cmpy), vi(_vi), N(_N), 
      parity(_parity) 
{
  setwidth2D();
}

template< typename T, typename CX, typename CY >
void d2sortvec<T,CX,CY>::setwidth2D()
{
  width = (int)(sqrt(N)); 
}

template< typename T, typename CX, typename CY >
void d2sortvec<T,CX,CY>::setwidth3D()
{
  width = (int)pow(N,2.0/3.0);
}


template< typename T, typename CX, typename CY >
void d2sortvec<T,CX,CY>::eval()
{
  //assert(width!=0);
  assert(N!=0);

  sort(&vi[0],&vi[0]+N,cmpx);

  uint k=parity;
  for (uint i=0; i<N; i+=width, ++k)
  {
    if (i+width<=N)
    {
      if (k%2)
        sort(vi+i,vi+i+width,mynot2(cmpy));
      else
        sort(vi+i,vi+i+width,cmpy);
    }
    else
    {
      if (k%2)
        sort(vi+i,vi+N,mynot2(cmpy));
      else
        sort(vi+i,vi+N,cmpy);
    }
  }
}





#endif