/****************************************************************************/
/* This overloaded class (a) organizes a matrix in order of 
	increasing values - uses the quicksort algorithm,
	(b) constructs an index table. 
	C.A. Bertulani - 03/09/2000    */
/****************************************************************************/

#include"butil.h"

class sort{
public:
	sort(unsigned long n, Number arr[]);  /* ordering */
	sort(unsigned long n, Number arr[], unsigned long indx[]); /* indexing */
	sort(unsigned long n, Number ra[], Number rb[], Number rc[]);
	sort(unsigned long n, Number a[], unsigned long indx[], unsigned long irank[]);
private:
	static const int M;
	static const int NSTACK;
	static const int FM;
	static const int FA;
	static const int FC;
};

const int sort::M = 7;
const int sort::NSTACK = 50;
/* constants for little random number generator  */
const int sort::FM = 7875;
const int sort::FA = 211;
const int sort::FC = 1663;
/**/

/* Constructor 'sort': orders a matrix  */
sort::sort(unsigned long n, Number arr[]){


	unsigned long l=1,jstack=0,j,ir,iq,i;
	int* istack;
//	int* istack = new int[NSTACK+1];
	istack=ivector(NSTACK);

	long int fx=0L;
	Number a; 

	ir=n;
	for (;;) {
		if (ir-l < M) { 
			/* sort by insertion for small subarray  */
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				for (i=j-1;arr[i]>a && i>=1;i--) arr[i+1]=arr[i];
				arr[i+1]=a;
			}
			if (jstack == 0) return;
			/**/
			/*  begin a new round of partitioning  */
			ir=istack[jstack--];   
			l=istack[jstack--];
		} else {
			i=l;
			j=ir;
			/*  generate a random element as a partitioning element  */
			fx=(fx*FA+FC) % FM;     
			iq=l+((ir-l+1)*fx)/FM;
			/**/
			a=arr[iq];
			arr[iq]=arr[l];
			for (;;) {
				/*  scan down to find element < a  */
				while (j > 0 && a < arr[j]) j--;
				if (j <= i) { arr[i]=a;	break; } /* pointers crossed */
				arr[i++]=arr[j];
				/*  scan up to find element > a  */
				while (a > arr[i] && i <= n) i++;
				if (j <= i) { arr[(i=j)]=a; break; } /* pointers crossed */
				arr[j--]=arr[i];
			}
			/* push pointers to larger subarray on stack, 
			   process smaller subarray immediately */
			if (ir-i >= i-l) {
				istack[++jstack]=i+1;
				istack[++jstack]=ir;
				ir=i-1;
			} else {
				istack[++jstack]=l;
				istack[++jstack]=i-1;
				l=i+1;
			}
			/* handle exceptions */
			try {if (jstack > NSTACK){ 
				throw "NSTACK too small in SORT";
			}
			}        
			catch (char* message) {
            cout << message;
			}
		}
	}
	delete [] istack;
}
/* end first constructor  */


/* Another overloaded constructor 'sort': indexes a set of records.  
   Indexes an array arr[1..n], i.e., outputs the array indx[1..n]
   such that arr[indx[j]] is in ascending order for j=1,2,..,N. 
   I.e., the integer initially numbered the smallest key ends up in 
   the number one position, and so on.
   The input quantities n and arr are not changed.  */

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


sort::sort(unsigned long n, Number arr[], unsigned long indx[])
{
	unsigned long i,indxt,ir=n,itemp,j,k,l=1;
	int jstack=0,*istack;
	Number a;
	istack=ivector(NSTACK);
	for (j=1;j<=n;j++) indx[j]=j;
	for (;;) {
		if (ir-l < M) {
			/* sort by insertion for smal subarray */
			for (j=l+1;j<=ir;j++) {
				indxt=indx[j];
				a=arr[indxt];
				for (i=j-1;i>=1;i--) {
					if (arr[indx[i]] <= a) break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if (jstack == 0) break;
			/**/
			/* begin a new round of partitioning */
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			/* choose median of left, center and right elements
			as partitioning element a. Also rearrange so that 
			a[l] <= a[l+1] <= a[ir]
			*/
			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])
			}
			/* initialize pointers */
			i=l+1;
			j=ir;
			/* partitioning element */
			indxt=indx[l+1];
			a=arr[indxt];
			/* scan up and down to find element > a and < a, respectively  */
			for (;;) {
				do i++; while (arr[indx[i]] < a);
				do j--; while (arr[indx[j]] > a);
				if (j < i) break;      /* pointers crossed */
				SWAP(indx[i],indx[j])    /* exchange elements */
			}
			indx[l+1]=indx[j];
			indx[j]=indxt;      /* insert partioning element */
			jstack += 2;
			/* handle exceptions */
			try {if (jstack > NSTACK){ 
				throw "NSTACK too small in SORT";
			}
			}        
			catch (char* message) {
            cout << message;
			}
			/* push pointers to larger subarray on stack, 
			   process smaller subarray immediately */
			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;
			}
		}
	}
	delete [] istack;
}
#undef SWAP

/* end second constructor  */


/* Another overloaded constructor 'sort': order a matriz
   with corresponding rearrangementes in two other.  */

sort::sort(unsigned long n, Number ra[], Number rb[], Number rc[])

/* sorts an array ra[1..n] into ascending numerical order while making
the corresponding rearrangements of the arryas rb[1..n] and rc[1..n].
An index table is constructed via the routine sort ( #, #, #).
*/
{
	unsigned long j,*iwksp;
	Number *wksp;
	iwksp = lvector(n);
	wksp = vector(n);
	sort(n, ra, iwksp);    /*  make an index table */
	for (j=1;j<=n;j++) wksp[j]=ra[j];     /*  save the array ra */
	for (j=1;j<=n;j++) ra[j]=wksp[iwksp[j]];  /* copy it back in rearranged order */
	for (j=1;j<=n;j++) wksp[j]=rb[j];         /* same for the other two   */
	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]];
    delete [] iwksp;
	delete [] wksp;
}

/* end third constructor  */


sort::sort(unsigned long n, Number a[], unsigned long ind[], unsigned long irank[])
{
	unsigned long j;

	sort(n, a, ind);    /*  make an index table */
	for (j=1;j<=n;j++) irank[ind[j]]=j;
}

#undef M
#undef NSTACK
#undef FM
#undef FA
#undef FC

/*    End class sort   */