/* 
   Copyright (C) 1999 E. H. Haley

   top.c: 
   
   top (mineminemine) calls NR functions in sel.c:

   sorti (sort w/ index table)
   sort2 (quicksort?)
   which both call indexx

   For that reason, unfortunately, top too is NRly enough that you
   need to send in arr-1, heap-1, index-1 to account for their taste.
   EWW!  this also means makes the indices they return are ALL off by
   one, forever!!  (All the way down.)  Therefore the last thing top.c
   does is decrement all of them.  "One kludge duct-taped to another."

*/

#include "top.h"
#include "sel.h"

void top(unsigned long m, unsigned long n, float arr[], float heap[],
	  unsigned long index[])

     /*Returns in heap[1..m] the largest m elements of the array
       arr[1..n], in descending absval order, and in index[] their
       original indices. The array arr is not altered. For e ciency,
       this routine should be used only when m<<n.*/

     /*not only do I have fi&ffi deligatureized, also shall I have an
       index array */
{ 
  void sorti(unsigned long n, float arr1[], unsigned long arr2[]); 
  void sort2(unsigned long n, float arr1[], unsigned long arr2[]); 
  void nrerror(char error_text[]); 
  unsigned long i,j,k; 
  float swap;

  for (i=1;i<=m;i++)
    {
      index[i]=i;
      heap[i]=arr[i];  
    }

  sorti(m,heap,index);
  //Create initial heap by overkill! We assume m<<n. 
  for (i=m+1;i<=n;i++)       //For each remaining element... 
    { 
      //Put it on the heap?
      if ((arr[i]>0?arr[i]:-arr[i]) > (heap[1]>0?heap[1]:-heap[1]))
	{ 
	  heap[1]=arr[i]; 
	  index[1]=i;
	  for (j=1;;) 	      //Sift down.
	    { 
	      k=j << 1; 

	      if (k > m) break; 
	      if (k != m && 
		  ((heap[k]>0?heap[k]:-heap[k])
		   > (heap[k+1]>0?heap[k+1]:-heap[k+1]))
		  ) k++; 
	      if (
		  (heap[j]>0?heap[j]:-heap[j]) 
		  <= (heap[k]>0?heap[k]:-heap[k])
		  ) break; 
	      
	      swap=heap[k]; 
	      heap[k]=heap[j]; 
	      heap[j]=swap;

	      swap=index[k];
	      index[k]=index[j];
	      index[j]=swap;
	      j=k; 
	    } 
	} 
    } 
  sort2(m,heap,index);
  for (i=1; i<=m; i++)
     index[i]--;
}

/*
#include <stdlib.h>
#include <stdio.h>
int main(void)
{
  unsigned long i;
  float arr[100], heap[10];
  unsigned long index[10];

  for(i=0; i<100; i++)
    arr[i]= (float)random()/((float)RAND_MAX+1.0)-0.5;
  top(10,100, arr-1, heap-1, index-1);
  for(i=0; i<100; i++)    printf("%ld\t%.3f\n",i,arr[i]);
  for(i=0; i<10; i++)
    printf("%ld\t%ld\t%.3f\n",i,index[i],arr[index[i]]);

  exit(0);
}
*/