# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "hilbert_curve.h"

/******************************************************************************/

void d2xy ( int m, int d, int *x, int *y )

/******************************************************************************/
/*
  Purpose:

    d2xy() converts a 1D Hilbert coordinate to a 2D Cartesian coordinate.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    05 December 2015

  Input:

    int M, the index of the Hilbert curve.
    The number of cells is N=2^M.
    0 < M.

    int D, the Hilbert coordinate of the cell.
    0 <= D < N * N.

  Output:

    int *X, *Y, the Cartesian coordinates of the cell.
    0 <= *X, *Y < N.
*/
{
  int n;
  int rx;
  int ry;
  int s;
  int t = d;

  n = i4_power ( 2, m );

  *x = 0;
  *y = 0;
  for ( s = 1; s < n; s = s * 2 )
  {
    rx = 1 & ( t / 2 );
    ry = 1 & ( t ^ rx );
    rot ( s, x, y, rx, ry );
    *x = *x + s * rx;
    *y = *y + s * ry;
    t = t / 4;
  }
  return;
}
/******************************************************************************/

int i4_power ( int i, int j )

/******************************************************************************/
/*
  Purpose:

    i4_power() returns the value of I^J.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 October 2007

  Author:

    John Burkardt

  Input:

    int I, J, the base and the power.  J should be nonnegative.

  Output:

    int I4_POWER, the value of I^J.
*/
{
  int k;
  int value;

  if ( j < 0 )
  {
    if ( i == 1 )
    {
      value = 1;
    }
    else if ( i == 0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "i4_power(): Fatal error!\n" );
      fprintf ( stderr, "  I^J requested, with I = 0 and J negative.\n" );
      exit ( 1 );
    }
    else
    {
      value = 0;
    }
  }
  else if ( j == 0 )
  {
    if ( i == 0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "i4_power(): Fatal error!\n" );
      fprintf ( stderr, "  I^J requested, with I = 0 and J = 0.\n" );
      exit ( 1 );
    }
    else
    {
      value = 1;
    }
  }
  else if ( j == 1 )
  {
    value = i;
  }
  else
  {
    value = 1;
    for ( k = 1; k <= j; k++ )
    {
      value = value * i;
    }
  }
  return value;
}
/******************************************************************************/

void rot ( int n, int *x, int *y, int rx, int ry ) 

/******************************************************************************/
/*
  Purpose:

    rot() rotates and flips a quadrant appropriately.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    05 December 2015

  Input:

    int N, the length of a side of the square.  N must be a power of 2.

    int *X, *Y, the coordinates of a point.

    int RX, RY, values of 0 or 1 which indicate
    whether reflection or flips should be carried out.

  Output:

    int *X, *Y, the transformed coordinates.
*/
{
  int t;

  if ( ry == 0 )
  {
/*
  Reflect.
*/
    if ( rx == 1 )
    {
      *x = n - 1 - *x;
      *y = n - 1 - *y;
    }
/*
  Flip.
*/
     t = *x;
    *x = *y;
    *y =  t;
  }
  return;
}
/******************************************************************************/

int xy2d ( int m, int x, int y )

/******************************************************************************/
/*
  Purpose:

    XY2D converts a 2D Cartesian coordinate to a 1D Hilbert coordinate.

  Discussion:

    It is assumed that a square has been divided into an NxN array of cells,
    where N is a power of 2.

    Cell (0,0) is in the lower left corner, and (N-1,N-1) in the upper 
    right corner.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    05 December 2015

  Input:

    int M, the index of the Hilbert curve.
    The number of cells is N=2^M.
    0 < M.

    int X, Y, the Cartesian coordinates of a cell.
    0 <= X, Y < N.

  Output:

    int XY2D, the Hilbert coordinate of the cell.
    0 <= D < N * N.
*/
{
  int d = 0;
  int n;
  int rx;
  int ry;
  int s;

  n = i4_power ( 2, m );

  for ( s = n / 2; s > 0; s = s / 2 )
  {
    rx = ( x & s ) > 0;
    ry = ( y & s ) > 0;
    d = d + s * s * ( ( 3 * rx ) ^ ry );
    rot ( s, &x, &y, rx, ry );
  }
  return d;
}