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

#define MAXN 100000

void main()
{

  double value[MAXN],checksum;
  int    n,seed,m,i;
  void   initialize_values(),hpsort();

/* User inputs. */

  n = MAXN + 1;
  while (n > MAXN) {
    printf("Enter number of values to sort: \n");
    scanf("%d", &n);
  }

  printf("Enter random number seed: \n");
  scanf("%d", &seed);

  m = n + 1;
  while (m > n) {
    printf("Enter which final value to print out: \n");
    scanf("%d", &m);
  }

/* Initialize the values. */

  initialize_values(n,value,seed);

/* Sort the values. 
   Note: Numerical Recipes convention is value(1:n), so must call
         hpsort with value-1
*/

  hpsort(n,value-1);

/* Compute and print output. */

  checksum = 0.0;
  for (i = 0; i < n; i++) checksum += (i+1)*value[i];
  checksum = checksum/n;

  printf("Sorted value %d = %g\n",m,value[m-1]);
  printf("Checksum = %g\n",checksum);

}


/* Initialize the values randomly. */

void initialize_values(n,value,seed)
  int    n;
  double value[];
  int    seed;
{

  int    i;
  double rand1();

/* generate random values
   random number generator is fully periodic, so will not create duplicates
*/

  for (i = 0; i < n; i++)
    value[i] = rand1(&seed);
}


/* Heap sort from Numerical Recipes. */

void hpsort(n,ra)
  int n;
  double *ra;
{
  int i,ir,j,l;
  double rra;
  
  if (n < 2) return;
  l=(n >> 1)+1;
  ir=n;
  for (;;) {
    if (l > 1) {
      rra=ra[--l];
    } else {
      rra=ra[ir];
      ra[ir]=ra[1];
      if (--ir == 1) {
	ra[1]=rra;
	break;
      }
    }
    i=l;
    j=l+l;
    while (j <= ir) {
      if (j < ir && ra[j] < ra[j+1]) j++;
      if (rra < ra[j]) {
	ra[i]=ra[j];
	i=j;
	j <<= 1;
      } else j=ir+1;
    }
    ra[i]=rra;
  }
}


/* Park/Miller RNG */
 
double rand1(iseed)
  int *iseed;
{

  double aa=16807.0;
  double mm=2147483647.0;
  double sseed;
  int    jseed;

  jseed = *iseed;
  sseed = jseed;
  jseed = aa*sseed/mm;
  sseed = aa*sseed - mm*jseed;
  jseed = sseed;
  
  *iseed = jseed;

  return(sseed/mm);

}