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

/* An implementation of heapsort as an index sorting algorithm for
   double precision arrays. */

/* Swap IX[a] and IX[b]. */
void swap(int* IX,
          const int a,
          const int b)
{
  int U = IX[a];
  IX[a] = IX[b];
  IX[b] = U;
}


/* Sift node I of V down through the tree. */
void sift_d(const double* V,
	    int* IX,
	    const int I,
	    const int n)
{
  /* Store the new location of node I in ii. */
  int ii = I;

  /* The lowest node of V with a child. */
  int IL = n/2 - 1;

  while (1) 
    {
      /* Indices of the children. */
      int j1 = 2*ii + 1;
      int j2 = (j1 + 1 >= n) ? j1 : j1 + 1;
      
      /* Index of the greater child. */
      int g = (V[IX[j1]] > V[IX[j2]]) ? j1 : j2;
      
      /* Push down the parent if the child is greater. */
      if (V[IX[ii]] < V[IX[g]]) 
        {
          swap(IX, ii, g);
          ii = g;

	  /* We've reached the bottom. */
          if (g > IL) break;
        }
      
      /* Otherwise we've pushed as far as we can go. */
      else break;
    }
}


/* Sift node I of V down through the tree.  This is the integer
   version. */
void sift_i(const int* V,
	    int* IX,
	    const int I,
	    const int n)
{
  int ii, j1, j2, g;

  /* Store the new location of node I in ii. */
  ii = I;

  /* The lowest node of V with a child. */
  int IL = n/2 - 1;

  while (1) 
    {
      /* Indices of the children. */
      j1 = 2*ii + 1;
      j2 = (j1 + 1 >= n) ? j1 : j1 + 1;
      
      /* Index of the greater child. */
      g = (V[IX[j1]] > V[IX[j2]]) ? j1 : j2;
      
      /* Push down the parent if the child is greater. */
      if (V[IX[ii]] < V[IX[g]]) 
        {
          swap(IX, ii, g);
          ii = g;
          if (g > IL) break;
        }
      
      /* Otherwise we've pushed as far as we can go. */
      else break;
    }
}


/* Index sort a double array using heapsort. */
void heap_index_sort_d(const double* Z,
		       int* IX,
		       const int n) 
{
  int i, J;

  /* Initialize IX to the indentity. */
  for (i=0; i<n; ++i) IX[i] = i;

  /* Make a heap. */
  for (i=n/2-1; i>=0; --i) sift_d(Z, IX, i, n);

  /* Sort the heap. */ 
  for (i=n-1; i>0; --i) 
    {
      /* Pull of the top node. */
      J = IX[0];
      
      /* Move the lowest node to the top, push it down to its place. */
      IX[0] = IX[i];
      sift_d(Z, IX, 0, i);

      /* Put the old top node into the sorted part of the array. */
      IX[i] = J;
    }
}


/* Index sort an integer array using heapsort. */
void heap_index_sort_i(const int* Z,
		       int* IX,
		       const int n) 
{
  int i, J;

  /* Initialize IX to the indentity. */
  for (i=0; i<n; ++i) IX[i] = i;

  /* Make a heap. */
  for (i=n/2-1; i>=0; --i) sift_i(Z, IX, i, n);

  /* Sort the heap. */ 
  for (i=n-1; i>0; --i) 
    {
      /* Pull of the top node. */
      J = IX[0];
      
      /* Move the lowest node to the top, push it down to its place. */
      IX[0] = IX[i];
      sift_i(Z, IX, 0, i);

      /* Put the old top node into the sorted part of the array. */
      IX[i] = J;
    }
}


/* Compute the ranks for the elements in Z. */
void rank(const double* Z,
	  int* R,
	  const int n)
{
  int* IX = (int*)malloc(n*sizeof(int));
  heap_index_sort_d(Z, IX, n);
  heap_index_sort_i(IX, R, n);
  free(IX);
}


int main(int argc, 
         char** argv)
{
  int L, k;
  double *V;
  int *IX, *R;

  srand48(7348348);

  if (argc != 2)
    {
      printf("Usage: %s <length>\n", argv[0]);
      exit(0);
    }

  /* Construct a random array. */
  L = atoi(argv[1]);
  V = (double*)malloc(L*sizeof(double));
  for (k=0; k<L; ++k) V[k] = drand48();

  /* Storage for the indices and the ranks. */
  IX = (int*)malloc(L*sizeof(int));
  R = (int*)malloc(L*sizeof(int));

  /* Do the sort. */
  heap_index_sort_d(V, IX, L);
  
  /* Get the ranks. */
  heap_index_sort_i(IX, R, L);
  
  /* Print the sorted array. */
  for (k=0; k<L; ++k) printf("%12.4f\t%4d\t%4d\n", V[k], IX[k], R[k]);

  return 0;
}