sort.c

地震波正演和显示模块

/*********************** self documentation **********************/
/*****************************************************************************
SORT - Functions to sort arrays of data or arrays of indices

hpsort		sort an array of floats by the heap sort method
qkisort		sort an array of indices i[] so that 
		a[i[0]] <= a[i[1]] <= ... <= a[i[n-1]]
		uses the "quick sort" method
qkifind		partially sort an array of indices i[] so that the 
		index i[m] has the value it would have if the entire
		array of indices were sorted such that 
		a[i[0]] <= a[i[1]] <= ... <= a[i[n-1]]
		uses the "quick sort" method
qksort		sort an array of floats such that a[0] <= a[1] <= ... <= a[n-1]
		uses the "quick sort" method
qkfind		partially sort an array of floats  so that the element a[m] has
		the value it would have if the entire array were sorted
		such that a[0] <= a[1] <= ... <= a[n-1]
		uses the "quick sort" method

******************************************************************************
Function Prototypes:
void hpsort (int n, float a[]);
void qkisort (int n, float a[], int i[]);
void qkifind (int m, int n, float a[], int i[]);
void qksort (int n, float a[]);
void qkfind (int m, int n, float a[]);

******************************************************************************
hpsort:
Input:
n		number of elements in a
a		array[n] to be sorted

Output:
a		array[n] sorted

******************************************************************************
qkisort:
Input:
n		number of elements in array a
a		array[n] elements
i		array[n] indices to be sorted

Output:
i		array[n] indices sorted

******************************************************************************
qkifind:
Input:
m		index of element to be found
n		number of elements in array a
a		array[n] elements
i		array[n] indices to be partially sorted

Output:
i		array[n] indices partially sorted sorted

******************************************************************************
qksort:
Input:
n		number of elements in array a
a		array[n] containing elements to be sorted

Output:
a		array[n] containing sorted elements

******************************************************************************
qkfind:
Input:
m		index of element to be found
n		number of elements in array a
a		array[n] to be partially sorted

Output:
a		array[n] partially sorted

*****************************************************************************/
/**************** end self doc ********************************/


void
hpsort (int n, float a[])
/*****************************************************************************
sort an array so that a[0] <= a[1] <= ... <= a[n-1]
******************************************************************************
Input:
n		number of elements in a
a		array[n] to be sorted

Output:
a		array[n] sorted

*****************************************************************************/
{
	int kroot,klast,kparent,kchild;
	float aparent;
	
	/* initialize indices of root node and last node in heap */
	kroot = n/2-1;
	klast = n-1;
	
	/* loop until array is sorted */
	while (klast>0) {
	
		/* if in the heap building phase */
		if (kroot>=0) {
		
			/* set index of parent to be added to the heap */
			kparent = kroot;

			/* set value of parent to be added to the heap */
			aparent = a[kroot--];
		
		/* else, the tree is a heap and in the sorting phase */
		} else {
			
			/* set index of parent at which to begin sifting */
			kparent = 0;
	
			/* set parent value to last value in heap */
			aparent = a[klast];
			
			/* copy top of heap to sorted elements at end */
			a[klast--] = a[0];
		}

		/* sift parent down until greater than both children */
		for (kchild=kparent*2+1; kchild<=klast; kchild=kparent*2+1) {

			/* get index of greater of two children */
			if (kchild<klast && a[kchild+1]>a[kchild]) kchild++;

			/* if greater child is greater than parent */
			if (a[kchild]>aparent) {
			
				/* promote the greater child */
				a[kparent] = a[kchild];

				/* demote the parent */
				kparent = kchild;
			
			/* else, if parent is greater than children, break */
			} else 
				break;
		}
	
		/* set value of parent (which may have been demoted) */
		a[kparent] = aparent;
	}
}



#define NSTACK 50	/* maximum sort length is 2^NSTACK */
#define NSMALL 7	/* size of array for which insertion sort is fast */
#define FM 7875		/* constants used to generate random pivots */
#define FA 211
#define FC 1663

static void
qkipart (float a[], int i[], int p, int q, int *j, int *k)
/*****************************************************************************
quicksort partition (FOR INTERNAL USE ONLY):
take the value x of a random element from the subarray a[p:q] of
a[0:n-1] and rearrange indices in the subarray i[p:q] in such a way
that there exist integers j and k with the following properties:
  p <= j < k <= q, provided that p < q
  a[i[l]] <= x,  for p <= l <= j
  a[i[l]] == x,  for j < l < k
  a[i[l]] >= x,  for k <= l <= q
note that this effectively partitions the subarray with bounds
[p:q] into lower and upper subarrays with bounds [p:j] and [k:q]
******************************************************************************
Input:
a		array[p:q]
i		array[p:q] of indices to be rearranged
p		lower bound of subarray; must be less than q
q		upper bound of subarray; must be greater then p

Output:
i		array[p:q] of indices rearranged
j		upper bound of lower output subarray
k		lower bound of upper output subarray

*****************************************************************************/
{
	int pivot,left,right,temp;
	float apivot;
	static long int seed=0L;
 
	/* choose random pivot element between p and q, inclusive */
	seed = (seed*FA+FC)%FM;
	pivot = p+(q-p)*(float)seed/(float)FM;
	if (pivot<p) pivot = p;
	if (pivot>q) pivot = q;
	apivot = a[i[pivot]];

	/* initialize left and right pointers and loop until break */
	for (left=p,right=q;;) {
		/*
		 * increment left pointer until either
		 * (1) an element greater than the pivot element is found, or
		 * (2) the upper bound of the input subarray is reached
		 */
		while (a[i[left]]<=apivot && left<q) left++;
 
		/*
		 * decrement right pointer until either
		 * (1) an element less than the pivot element is found, or
		 * (2) the lower bound of the input subarray is reached
		 */
		while (a[i[right]]>=apivot && right>p) right--;
 
		/* if left pointer is still to the left of right pointer */
		if (left<right) {

			/* exchange left and right indices */
			temp = i[left];
			i[left++] = i[right];
			i[right--] = temp;
		} 
		/* else, if pointers are equal or have crossed, break */
		else break;
	}
	/* if left pointer has not crossed pivot */
	if (left<pivot) {

		/* exchange indices at left and pivot */
		temp = i[left];
		i[left++] = i[pivot];
		i[pivot] = temp;
	}
	/* else, if right pointer has not crossed pivot */
	else if (pivot<right) {

		/* exchange indices at pivot and right */
		temp = i[right];
		i[right--] = i[pivot];
		i[pivot] = temp;
	}
	/* left and right pointers have now crossed; set output bounds */
	*j = right;
	*k = left;
}

static void
qkiinss (float a[], int i[], int p, int q)
/*****************************************************************************
quicksort insertion sort (FOR INTERNAL USE ONLY):
Sort a subarray of indices bounded by p and q so that
a[i[p]] <= a[i[p+1]] <= ... <= a[i[q]]
******************************************************************************
Input:
a		subarray[p:q] containing elements
i		subarray[p:q] containing indices to be sorted
p		lower bound of subarray; must be less than q
q		upper bound of subarray; must be greater then p

Output:
i		subarray[p:q] of indices sorted

*****************************************************************************/
{
	int j,k,ij;
	float aij;

	for (j=p+1; j<=q; j++) {
		for (ij=i[j],aij=a[ij],k=j; k>p && a[i[k-1]]>aij; k--)
			i[k] = i[k-1];
		i[k] = ij;
	}
}

void
qkisort (int n, float a[], int i[])
/*****************************************************************************
Sort an array of indices i[] so that 
a[i[0]] <= a[i[1]] <= ... <= a[i[n-1]]
******************************************************************************
Input:
n		number of elements in array a
a		array[n] elements
i		array[n] indices to be sorted

Output:
i		array[n] indices sorted

*****************************************************************************/
{
	int pstack[NSTACK],qstack[NSTACK],j,k,p,q,top=0;

	/* initialize subarray lower and upper bounds to entire array */