/*
  File name: fsort.c
  Created by: Ljubomir Buturovic
  Created: 05/02/2002
  Purpose: floating-point sorting functions based on code by
  R. Sedgewick.
*/

#include <float.h>
#include <stdlib.h>
#include "pcp.h"
#include "lau.h"
#include "sort.h"

static char rcsid[] = "$Id: fsort.c,v 1.12 2005/06/25 08:35:00 ljubomir Exp $";

/*
  In-place float sorting function based on algorithm by
  R. Sedgewick. In my experience, it is about 40% faster than the C
  library function qsort(). The sedgesort() function which implements
  the algorithm is downloaded from
  http://www.yendor.com/programming/sort/. ljb, 10/14/2001

  BIG NOTE: 'vector' needs to have len+1 elements, instead of the
  usual len. The first len elements are the numbers to be sorted. The
  last element is used internally by sedgesort() function (which is
  called by ssort()) as a 'sentinel'. This is obviously a restriction
  on part on finsort(); it means that the user must allocate one more
  element for each vector to be sorted. If this is not acceptable, use
  fsort(), which creates a new vector with the extra element.
*/
void finsort(float *vector, int len)
{
  if ((vector != (float *) 0) && (len > 0))
    {
      vector[len] = FLT_MAX;
      sedgesort(vector, len);
    }
}

/*
  Floating-point sorting function. The function allocates space for
  the sorted vector. The input 'vector' is not modifed.

  Return (float *) in case of error and set errno.
*/
float *fsort(float *vector, int len)
{
  float *vec;

  vec = malloc((len+1)*sizeof(float));
  if (vec != (float *) 0)
    {
      fvec_copy(vec, vector, len);
      finsort(vec, len);
    }
  return vec;
}

/*
  Floating-point sorting function. hsort() uses modified version of
  af_heapsort() to sort 'm' smallest values in a vector. This makes
  sense if 'k' is considerably smaller than the length of the vector;
  otherwise, sedgesort() is faster.
*/
void hsort(float *vector, int *len, int *m)
{
  int i;
  int length;
  int n;
  int limit;
  float temp;

  length = *len;
  n = *m;
  for (i = 0; i < length; i++)
    vector[i] = -vector[i];
  af_heapsort(vector, length, n);
  limit = length/2;
  if (n < limit)
    limit = n;
  for (i = 0; i < limit; i++)
    {
      temp = vector[i];
      vector[i] = vector[length-1-i];
      vector[length-1-i] = temp;
    }
  for (i = 0; i < n; i++)
    vector[i] = -vector[i];
}

/*
 |  void  af_heapsort (array, len)
 |  KEY_T  array[];
 |  int    len;
 |
 |  Abstract:	Sort array[0..len-1] into increasing order.
 |
 |  Method:	Heap-sort (ala Jon Bentley)
 |
 |  Advantages:
 |
 |	1. Stable sorting method (doesn't spoil relative order between
 |	   equal-key elements)
 |	2. Worst case complexity is O (n log(n))
 |	3. Method of choice when only the smallest n elemens are wanted
 |	   and not a full ordering (use sift n times).
 */

void  af_heapsort (KEY_T array[], int len, int n)
{
    register int	i;
    register KEY_T	temp, *sa = array - 1;
    int limit;
    /*
     |  'sa[]' is 'array[]' "shifted" one position to the left
     |  i.e.  sa[i] == array[i - 1]  (in particular:   sa[1] == array[0])
     |  'sa' has "Pascalish" indices and is thus more convenient
     |  for heap-sorting. An index i obeys the law:
     |	left_child (i) == 2i,  right_child (i) == 2i + 1
     |  Note: we don't use a[-1] which is outside our range, we just
     |  fake an array that starts one address lower, so we can refer
     |  to its first element as a[1] (rather than as a[0])
     */

    /* first step: make 'sa[]' a heap using siftdown len/2 times */
    for (i = len / 2; i >= 1; i--)
	siftdown (sa, i, len);

    /* heapify n times, to reach partial order */
    limit = 2;
    if (len-n+1 > 2)
      limit = len-n+1;
    for (i = len; i >= limit; i--) 
      {
	SWAP(sa[1], sa[i]);
	siftdown (sa, 1, i - 1);
      }
}


/*
 |  void  siftdown (array, lower, upper)
 |  KEY_T  array[];
 |  int    lower, upper;
 |
 |  Abstract:
 |	Assume array[] is a heap except array[lower] is not in order.
 |	Sift array[lower] down the heap as long as it is bigger than
 |	its two children nodes. Preserving the heap condition.
 |
 |	Optimization: rather than swapping array[i] with the greater
 |	of its children within the loop, we save array[lower]
 |	before the loop, shift the inner values within the loop, and
 |	finally plug the saved value into its final location.
 */
void  siftdown (KEY_T array[], int lower, int upper)
{
    register int	i = lower, c = lower;
    register int	lastindex = upper/2;
    register KEY_T	temp;

    temp = array[i];
    while (c <= lastindex) {
	c = 2 * i;		/* c = 2i is the left-child of i */
	if (c + 1 <= upper  &&  GT(array[c + 1], array[c]))
	    c++;		/* c is the greatest child of i */
	if (GE(temp, array[c])) break;
	array[i] = array[c];
	i = c;
    }
    array[i] = temp;
}