#include "/misc/radon2/xiaoxiak/thesis/header/nr.h"
void shell(unsigned long n, double a[])
{
/* Sorts an array a[] into ascending numerical order by Shell's method */
   unsigned long i,j,inc;
   double v;
   inc=1;
   do {
      inc *= 3;
      inc++;
   } while (inc <= n);
   do { 
      inc /= 3;
      for (i=inc+1;i<=n;i++) { 
         v=a[i];
         j=i;
         while (a[j-inc] > v) { 
            a[j]=a[j-inc];
            j -= inc;
	    if (j <= inc) break;
         }
         a[j]=v;
      }
   } while (inc > 1);
}

#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nr.h"
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50
                                
void sort(unsigned long n, double arr[])
{
/* quicksort for 1 vector */
        unsigned long i,ir=n,j,k,l=1,*istack;   
        int jstack=0;
        double a,temp;
                            
        istack=lvector(1,NSTACK);
        for (;;) {      
                if (ir-l < M) {
                        for (j=l+1;j<=ir;j++) {
                                a=arr[j];
                                for (i=j-1;i>=l;i--) {
                                        if (arr[i] <= a) break;
                                        arr[i+1]=arr[i];
                                }
                        }
                        if (jstack == 0) break;
                        ir=istack[jstack--];
                        l=istack[jstack--];
                } else { 
                        k=(l+ir) >> 1;
                        SWAP(arr[k],arr[l+1])
                        if (arr[l] > arr[ir]) {
                                SWAP(arr[l],arr[ir])
                        }
                        if (arr[l+1] > arr[ir]) {
                                SWAP(arr[l+1],arr[ir])
                        }
                        if (arr[l] > arr[l+1]) {
                                SWAP(arr[l],arr[l+1])
                        }
                        i=l+1;
                        j=ir;
                        a=arr[l+1];
                        for (;;) {
                                do i++; while (arr[i] < a);
                                do j--; while (arr[j] > a);
                                if (j < i) break;
                                SWAP(arr[i],arr[j]);
                        }
                        arr[l+1]=arr[j];
                        arr[j]=a;
                        jstack += 2;
                        if (jstack > NSTACK) nrerror("NSTACK too small in sort.");
                        if (ir-i+1 >= j-l) {
                                istack[jstack]=ir;
                                istack[jstack-1]=i;
                                ir=j-1;
                        } else {
                                istack[jstack]=j-1;
                                istack[jstack-1]=l;
                                l=i;
                        }
                }
        }
        free_lvector(istack,1,NSTACK);
}
#undef M
#undef NSTACK
#undef SWAP
#undef NRANSI


void hpsort(unsigned long n, double ra[])
{
/* Heapsort */
	unsigned long 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 break;
		}
		ra[i]=rra;
	}
}

#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define SWAP(a,b) itemp=(a);(a)=(b);(b)=itemp;
#define M 7
#define NSTACK 5000000
void indexx(unsigned long n, double arr[],unsigned long indx[])
{
/* construct index for array */
	unsigned long i,indxt,ir=n,itemp,j,k,l=1;
	int jstack=0,*istack;
	double a;
	istack=ivector(1,NSTACK);
	for (j=1;j<=n;j++) indx[j]=j;
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				indxt=indx[j];
				a=arr[indxt];
				for (i=j-1;i>=l;i--) {
					if (arr[indx[i]] <= a) break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP(indx[k],indx[l+1]);
			if (arr[indx[l]] > arr[indx[ir]]) {
				SWAP(indx[l],indx[ir])
			}
			if (arr[indx[l+1]] > arr[indx[ir]]) {
				SWAP(indx[l+1],indx[ir])
			}
			if (arr[indx[l]] > arr[indx[l+1]]) {
				SWAP(indx[l],indx[l+1])
			}
			i=l+1;
			j=ir;
			indxt=indx[l+1];
			a=arr[indxt];
			for (;;) {
				do i++; while (arr[indx[i]] < a);
				do j--; while (arr[indx[j]] > a);
				if (j < i) break;
				SWAP(indx[i],indx[j])
			}
			indx[l+1]=indx[j];
			indx[j]=indxt;
			jstack += 2;
			if (jstack > NSTACK) nrerror("NSTACK too small in indexx.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	free_ivector(istack,1,NSTACK);
}
#undef M
#undef NSTACK
#undef SWAP
#undef NRANSI


#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"

void sort3(unsigned long n, double ra[], double rb[], double rc[])
{
/* sort, use an index to sort 3 or more arrays */
	void indexx(unsigned long n, double arr[], unsigned long indx[]);
	unsigned long j,*iwksp;
	double *wksp;

	iwksp=lvector(1,n);
	wksp=dvector(1,n);
	indexx(n,ra,iwksp);
	for (j=1;j<=n;j++) wksp[j]=ra[j];
	for (j=1;j<=n;j++) ra[j]=wksp[iwksp[j]];
	for (j=1;j<=n;j++) wksp[j]=rb[j];
	for (j=1;j<=n;j++) rb[j]=wksp[iwksp[j]];
	for (j=1;j<=n;j++) wksp[j]=rc[j];
	for (j=1;j<=n;j++) rc[j]=wksp[iwksp[j]];
	free_dvector(wksp,1,n);
	free_lvector(iwksp,1,n);
}
#undef NRANSI


/* construct a rank table for an array */
void rank(unsigned long n, unsigned long indx[], unsigned long irank[])
{
	unsigned long j;

	for (j=1;j<=n;j++) irank[indx[j]]=j;
}


#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;

double select(unsigned long k, unsigned long n, double arr[])
{
/* find the kth smallest value in an array */
	unsigned long i,ir,j,l,mid;
	double a,temp;

	l=1;
	ir=n;
	for (;;) {
		if (ir <= l+1) {
			if (ir == l+1 && arr[ir] < arr[l]) {
				SWAP(arr[l],arr[ir])
			}
			return arr[k];
		} else {
			mid=(l+ir) >> 1;
			SWAP(arr[mid],arr[l+1])
			if (arr[l] > arr[ir]) {
				SWAP(arr[l],arr[ir])
			}
			if (arr[l+1] > arr[ir]) {
				SWAP(arr[l+1],arr[ir])
			}
			if (arr[l] > arr[l+1]) {
				SWAP(arr[l],arr[l+1])
			}
			i=l+1;
			j=ir;
			a=arr[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP(arr[i],arr[j])
			}
			arr[l+1]=arr[j];
			arr[j]=a;
			if (j >= k) ir=j-1;
			if (j <= k) l=i;
		}
	}
}
#undef SWAP


#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define M 64
#define BIG 1.0e30
#define FREEALL free_dvector(sel,1,M+2);free_lvector(isel,1,M+2);

double selip(unsigned long k, unsigned long n, double arr[])
{
/* finds the kth smallest value, without altering an array */
	void shell(unsigned long n, double a[]);
	unsigned long i,j,jl,jm,ju,kk,mm,nlo,nxtmm,*isel;
	double ahi,alo,sum,*sel;

	if (k < 1 || k > n || n <= 0) nrerror("bad input to selip");
	isel=lvector(1,M+2);
	sel=dvector(1,M+2);
	kk=k;
	ahi=BIG;
	alo = -BIG;
	for (;;) {
		mm=nlo=0;
		sum=0.0;
		nxtmm=M+1;
		for (i=1;i<=n;i++) {
			if (arr[i] >= alo && arr[i] <= ahi) {
				mm++;
				if (arr[i] == alo) nlo++;
				if (mm <= M) sel[mm]=arr[i];
				else if (mm == nxtmm) {
					nxtmm=mm+mm/M;
					sel[1 + ((i+mm+kk) % M)]=arr[i];
				}
				sum += arr[i];
			}
		}
		if (kk <= nlo) {
			FREEALL
			return alo;
		}
		else if (mm <= M) {
			shell(mm,sel);
			ahi = sel[kk];
			FREEALL
			return ahi;
		}
		sel[M+1]=sum/mm;
		shell(M+1,sel);
		sel[M+2]=ahi;
		for (j=1;j<=M+2;j++) isel[j]=0;

		for (i=1;i<=n;i++) {
			if (arr[i] >= alo && arr[i] <= ahi) {
				jl=0;
				ju=M+2;
				while (ju-jl > 1) {
					jm=(ju+jl)/2;
					if (arr[i] >= sel[jm]) jl=jm;
					else ju=jm;
				}
				isel[ju]++;
			}
		}
		j=1;
		while (kk > isel[j]) {
			alo=sel[j];
			kk -= isel[j++];
		}
		ahi=sel[j];
	}
}
#undef M
#undef BIG
#undef FREEALL
#undef NRANSI


void hpsel(unsigned long m, unsigned long n, double arr[], double heap[])
{
/* find M largest values, without altering an array */
	void sort(unsigned long n, double arr[]);
	void nrerror(char error_text[]);
	unsigned long i,j,k;
	double swap;

	if (m > n/2 || m < 1) nrerror("probable misuse of hpsel");
	for (i=1;i<=m;i++) heap[i]=arr[i];
	sort(m,heap);
	for (i=m+1;i<=n;i++) {
		if (arr[i] > heap[1]) {
			heap[1]=arr[i];
			for (j=1;;) {
				k=j << 1;
				if (k > m) break;
				if (k != m && heap[k] > heap[k+1]) k++;
				if (heap[j] <= heap[k]) break;
				swap=heap[k];
				heap[k]=heap[j];
				heap[j]=swap;
				j=k;
			}
		}
	}
}