#ifndef SNAKESORT_H 
#define SNAKESORT_H

#include <cmath>
using namespace std;

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


/*!
\brief Methods to order objects with the property that
       in a linear sense they are close together, in two
       and three dimensions.

Such orderings are important in incremental algorithms with
 a state where ording the points to be near each other
 gets an algorithmic speed up (constant access time).

There is an inherent instability in this sort as for large
 data sets the jump or change in a coordinate can be
 extremely large for a finite number of points. 
 Contrast/compare this with a Hilbert curve.
*/
class snakesort
{
public:

  /** Sort in blocks of width. */
  template< typename T, typename CMP >
  static void d1
  (
    T* vi,
    uintc N,
    uintc width,
    CMP compare,
    uintc parity=0
  );

  /** Sort on x-axis and in blocks on y-axis. */
  template< typename T, typename CMP >
  static void d2
  (
    T* vi,
    uintc N,
    CMP compare
  );

  /** Sort in x-axis, and in blocks on y-axis, and inside
      these blocks sort in smaller blocks on the z-axis. */
  template< typename T, typename CMP >
  static void d3
  (
    T* vi,
    uintc N,
    CMP compare
  );

};



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


/** Sort in blocks of width. */
template< typename T, typename CMP >
void snakesort::d1
(
  T* vi,
  uintc N,
  uintc width,
  CMP compare,
  uintc parity
)
{
  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(compare));
      else
        sort(vi+i,vi+i+width,compare);
    }
    else
    {
      if (k%2)
        sort(vi+i,vi+N,mynot2(compare));
      else
        sort(vi+i,vi+N,compare);
    }
  }
}

template< typename T, typename CMP >
void snakesort::d2
(
  T* vi,
  uintc N,
  CMP compare
)
{
  uint width=(int)sqrt(N);
  sort(vi,vi+N,mycomparexobject<CMP>());
  d1(vi,N,width,mycompareyobject<CMP>());
}

template< typename T, typename CMP >
void snakesort::d3
(
  T* vi,
  uintc N,
  CMP compare
)
{
  uint width=(int)pow(N,2.0/3.0);
  uint width2=(int)pow(N,1.0/3.0);
  sort(vi,vi+N,mycomparexobject<CMP>());
  d1(vi,N,width,mycompareyobject<CMP>());
  uint parity=0;
  mycomparezobject<CMP> compare2;
  for (uint i=0; i<N; i+=width)
  {
    if (i+width<N)
      d1(vi+i,width,width2,compare2,parity);
    else
      d1(vi+i,N-i,width2,compare2,parity);
    ++parity;
  }
}


#endif