#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#define N 10000 
#define SAMPLES 6 

long int counter;   /* this is a "global" counter  */
                    /* to be used in instruction   */
		    /* counting of the program     */

/********************************************************************/
/*                                                                  */
/*           Quicksort of an Array	                            */
/*                                                                  */
/********************************************************************/

void quicksort( int Q[], long int size ) {
    /* sort array Q into decreasing order ---  */
    /* max is at Q[0], min at last element.    */
    long int i,j;
    int p,x;

    if (size < 2) return;  /* no work to sort this */

    p = Q[0];  /* p is pivot value */
    
    i = 0; j = (size - 1);
    while (i <= j) 
	/* split Q into two parts,     */
        /* left part is > p,           */
        /* right part is <= p.         */
        if (Q[i]>p) 
           i++;
        else 
	   { x = Q[i]; Q[i] = Q[j]; Q[j] = x; j--; };

    /* As a result of the split, there are three */
    /* possibilities:  left part is emtpy,       */
    /*                 right part is empty,      */
    /*                 neither is empty.         */
    /* Need to insure, before recursion, that    */
    /* recursion is on smaller-size array than Q */
    
    if (i==0) 
       { /* left part is empty; this can only be */
         /* possible if p is maximum element of  */
         /* array Q --- meaning Q[0] was the max */
         /* which got swapped to Q[size-1].  So  */
         /* swap it back to Q[0] and recurse.    */
	 x = Q[0]; Q[0] = Q[size-1]; Q[size-1] = x;
         quicksort(Q+1,size-1);
         return;
       };

     /* if ((size-i)==0) */
         /* right part is empty; this implies    */
         /* that p is greater than all elements  */
         /* of Q --- which is impossible, since  */
         /* p is just a value from the array Q!  */

     /* both left and right parts are non-emtpy  */
     quicksort(Q,i);       /* sort left part  */
     quicksort(Q+i,size-i);  /* sort right part */

     return;
   
     };

int quicksel( int Q[], long int size, long int t ) {
    /* find the t-th largest in array Q  */
    /* by using quicksort and indexing;  */
    /* this is an O(n lg n) method.      */

    quicksort(Q,size);
    return( Q[t] );
    }

/********************************************************************/
/*                                                                  */
/*           Heap Sort of an Array	                            */
/*                                                                  */
/********************************************************************/

long int heapsize;  /* this is a "global" size of heap variable */

long int parent( long int p ) {
    /* return p's parent or -1 if p has no parent in heap */
    if (p<1) return -1; else return ((p-1)/2); };

long int childone( long int p ) {
    /* return p's left child or -1 if p has no left child */
    long int r;
    r = (2*p)+1;
    if (r>=heapsize) return -1; else return r; };

long int childtwo( long int p ) {
    /* return p's right child or -1 if p has no right child */
    long int r;
    r = (2*p)+2;
    if (r>=heapsize) return -1; else return r; };

void heapify( int H[], long int x ) {
    /* H is a heap except for position x */    
    long int r;
    if (x<=0) return;  /* but root is always OK */

    if (H[x] > H[parent(x)]) return;  /* it is heapified! */ 
    /* swap H[x] and H[parent(x)] */
       r = H[parent(x)];
       H[parent(x)] = H[x];
       H[x] = r;
    /* now heapify above */
    heapify( H, parent(x) );
    }; 
  
long int pullx( int H[], long int x ) {
    /* Eliminate x from the heap, and return   */
    /* the location of the "hole" that remains */ 
    /* by recursion to malong intain heap property  */

    /* if x is a leaf, there's no recursion */
    if (childone(x)<0 & childtwo(x)<0) return x;

    /* if x has only one child, it must be the */ 
    /* right child only - so promote right it  */
    if (childtwo(x)<0)
       { H[x] = H[childone(x)];
         return pullx( H, childone(x) ); };

    /* if x has two children, promote the one  */
    /* that is the smaller child               */
    if (H[childone(x)]<H[childtwo(x)])
       { H[x] = H[childone(x)];
         return pullx( H, childone(x) ); }
    else
       { H[x] = H[childtwo(x)];
         return pullx( H, childtwo(x) ); }
    };

int extract( int H[] ) {
    /* extract the root from heap H and reduce */
    /* the size of heap;  return extracted val */
    int r;
    long int p; 
    r = H[0]; 
    p = pullx( H, 0 );  

    /* now p is the "hole", which should be at */
    /* the last position instead -- so swap it */
    /* with the last position and re-heapify   */
    H[p] = H[heapsize - 1];
    heapsize = heapsize - 1;
    heapify( H, p ); 

    return r;
    };

void heapsort( int H[], long int size ) {
    /* sort array H into decreasing order ---  */
    /* max is at H[0], min at last element.    */
    long int i;

    /* phase one:  put everything into heap    */
    heapsize = 0;
    for (i=0; i<size; i++) 
       { heapsize++;   
         heapify( H, i); }; 

    /* phase two:  move min to end of H[..]    */
    for (i=(size-1); i>0; i--) 
       H[i] = extract( H );

    };

int heapsel( int A[], long int size, long int t ) {
    /* find the t-th largest in array A  */
    /* by using heapsort and indexing;   */
    /* this is an O(n lg n) method.      */

    heapsort(A,size);
    return( A[t] );
    }

/********************************************************************/
/*                                                                  */
/*          Find Kth Largest of Array 	                            */
/*                                                                  */
/********************************************************************/

void fivesort( int* a, int* b, int* c, int* d, int* e ) {
    /* sort five numbers (given by pointers)    */
    int i, sorted, x;
    int* A[5];

    /* copy input to array */
    A[0] = a; A[1] = b; A[2] = c; A[3] = d; A[4] = e;

    /* now just bubble sort array A */ 
    sorted = 0;
    while (sorted == 0)
      { sorted = 1;
        for (i=0; i<4; i++) 
           if (*A[i]>*A[i+1]) 
              { x = *A[i]; *A[i] = *A[i+1]; *A[i+1] = x;
                sorted = 0; };
      };
    };

void sortcols( int A[], long int size ) {
    /* sort five numbers, then put middle in    */
    /* the first position [it's the median]     */
    /* consider A to be split into five rows of */
    /* equal size (ignore remainder); then sort */
    /* each column of the five rows.            */
    long int i,j;

    /* compute length of a row  */
    j = size / 5;

    /* loop to sort each "column" of the rows   */
    for (i=0; i<j; i++)
       fivesort(A+(0*j)+i,A+(1*j)+i,A+(2*j)+i,A+(3*j)+i,A+(4*j)+i);
    };

int select( int A[], long int size, long int t );

int getpivot( int A[], long int size ) {
    /* find a good pivot value for finding the  */
    /* k-th largest using the recursive method  */
    long int i,j;

    /* compute length of a row */
    j = size / 5;
    if (j==0) return A[0];  /* if array small, exit */
    
    /* first split into five rows and sort cols */
    sortcols(A,size);

    /* return, as pivot, median of middle row */
    return select(A+2*j,j,j/2); 
    };

int select( int A[], long int size, long int t ) {
    /* find the t-th largest in an array */
    /* by the recursive O(n) method      */

    int p,x,sorted;
    long int i,j;

    /* first, eliminate nonsense cases */
    if (t<0) return A[0];
    if (t>(size-1)) return A[0];

    /* if there are 14 or fewer elements      */
    /* then avoid recursion -- just sort      */
    /* the elements and index for the answer  */
    /* [ the algorithm FAILS for small sizes  */
    /*   because 14/5 = 2 and the median of   */
    /*   rows of length 2 won't split array   */
    /*   into smaller units in all cases ]    */
    if (size<15) {
       sorted = 0;  /* use bubblesort to sort */
                    /* array in decreasing    */
                    /* order -- then index    */
       while (sorted == 0) {
	  sorted = 1;
          for (i=0; i<(size-1); i++) 
	     if (A[i]<A[i+1]) {
		x = A[i]; A[i] = A[i+1]; A[i+1] = x;
                sorted = 0;
                };
          }; 
       /* now just index to return value */
       return (A[t]);
       };

    /* find a good pivot value */
    p = getpivot( A, size ); 

    /* split A into "< p" and "> p" */
    i = 0; j = (size - 1);
    while (i<=j)  
       if (A[i]>p)
         { x = A[j]; A[j] = A[i]; A[i] = x;  j--; }
       else i++;  

    /* determine size of left and right parts */
    j = size - i;  /* j is # in right part    */
                   /* i is # in left part     */
                   /* A[j] is first in right  */ 

    /* determine whether desired value is in  */
    /* the right or the left part, and then   */
    /* use recursion to find that value       */
    /* [t=0 => max elem, t=(size-1) => min]   */    
    
    if (t<=(j-1)) 
       /* desired value in right part */
       return select(A+i,j,t);
    else
       /* desired value in left part  */
       return select(A,i,t-j);

    };

main()  {
   /* Declare variables */
   int A[N];
   int S;
   int median;
   long int n, i, j;

   /* Print starting message */
   printf("* * * Running %d Tests on N = %d * * *\n",SAMPLES,N);

   for (S=0; S<SAMPLES; S++) {

      /* Test quicksort method of median find   */
      /* This method is O(n lg n)               */

      /* Initialize A[] to random values        */
      /* Reset counter, test quicksort method   */
          counter = 0;
          for (i=0; i<N; i++)  A[i] = rand()%N;
          median = quicksel(A,N,N/2);

      /* Test heapsort method of median finding */
      /* This method is O(n lg n)               */

      /* Initialize A[] to random values        */
      /* Reset counter, test heapsort method    */
          counter = 0;
          for (i=0; i<N; i++)  A[i] = rand()%N;
          median = heapsel(A,N,N/2);

      /* Test recursive select method of median */
      /* finding.  This method is O(n)          */
   
      /* Initialize A[] to random values        */
      /* Reset counter, test selection method   */
          counter = 0;
          for (i=0; i<N; i++)  A[i] = rand()%N;
          median = select(A,N,N/2); 

      };

   }