1 #include "mrilib.h"
2
3 /*------------------------------------------------------------------------------*/
4 /*------------- insertion sort : sort an array of double in-place --------------*/
5
isort_double(int n,double * ar)6 static void isort_double( int n , double * ar )
7 {
8 register int j , p ; /* array indices */
9 register double temp ; /* a[j] holding place */
10 register double * a = ar ;
11
12 if( n < 2 || ar == NULL ) return ;
13
14 for( j=1 ; j < n ; j++ ){
15
16 if( a[j] < a[j-1] ){ /* out of order */
17 p = j ;
18 temp = a[j] ;
19 do{
20 a[p] = a[p-1] ; /* at this point, a[p-1] > temp, so move it up */
21 p-- ;
22 } while( p > 0 && temp < a[p-1] ) ;
23 a[p] = temp ; /* finally, put temp in its place */
24 }
25 }
26 return ;
27 }
28
29 #undef QS_CUTOFF
30 #undef QS_SWAP
31 #undef QS_SWAPI
32 #define QS_CUTOFF 20 /* cutoff to switch from qsort to isort */
33 #define QS_SWAP(x,y) (temp=(x), (x)=(y),(y)=temp )
34 #define QS_SWAPI(x,y) (itemp=(x),(x)=(y),(y)=itemp)
35 #ifndef QS_STACK
36 # define QS_STACK 9999
37 #endif
38
39 static int stack[QS_STACK] ; /* stack for qsort "recursion" */
40
41 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
42
qsrec_double(int n,double * ar,int cutoff)43 static void qsrec_double( int n , double * ar , int cutoff )
44 {
45 register int i , j ; /* scanning indices */
46 register double temp , pivot ; /* holding places */
47 register double * a = ar ;
48
49 int left , right , mst , itemp,nnew ;
50
51 /* return if too short (insertion sort will clean up) */
52
53 if( cutoff < 3 ) cutoff = 3 ;
54 if( n < cutoff || ar == NULL ) return ;
55
56 /* initialize stack to start with whole array */
57
58 stack[0] = 0 ;
59 stack[1] = n-1 ;
60 mst = 2 ;
61
62 /* loop while the stack is nonempty */
63
64 while( mst > 0 ){
65 right = stack[--mst] ; /* work on subarray from left -> right */
66 left = stack[--mst] ;
67
68 i = ( left + right ) / 2 ; /* middle of subarray */
69
70 /*----- sort the left, middle, and right a[]'s -----*/
71
72 if( a[left] > a[i] ) QS_SWAP( a[left] , a[i] ) ;
73 if( a[left] > a[right] ) QS_SWAP( a[left] , a[right] ) ;
74 if( a[i] > a[right] ) QS_SWAP( a[right] , a[i] ) ;
75
76 pivot = a[i] ; /* a[i] is the median-of-3 pivot! */
77 a[i] = a[right] ;
78
79 i = left ; j = right ; /* initialize scanning */
80
81 /*----- partition: move elements bigger than pivot up and elements
82 smaller than pivot down, scanning in from ends -----*/
83
84 do{
85 for( ; a[++i] < pivot ; ) ; /* scan i up, until a[i] >= pivot */
86 for( ; a[--j] > pivot ; ) ; /* scan j down, until a[j] <= pivot */
87
88 if( j <= i ) break ; /* if j meets i, quit */
89
90 QS_SWAP( a[i] , a[j] ) ;
91 } while( 1 ) ;
92
93 /*----- at this point, the array is partitioned -----*/
94
95 a[right] = a[i] ; a[i] = pivot ; /* restore the pivot */
96
97 /*----- signal the subarrays that need further work -----*/
98
99 nnew = 0 ;
100 if( (i-left) > cutoff ){ stack[mst++] = left; stack[mst++] = i-1 ; nnew++; }
101 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right; nnew++; }
102
103 /* if just added two subarrays to stack, make sure shorter one comes first */
104
105 if( nnew == 2 && stack[mst-3] - stack[mst-4] > stack[mst-1] - stack[mst-2] ){
106 QS_SWAPI( stack[mst-4] , stack[mst-2] ) ;
107 QS_SWAPI( stack[mst-3] , stack[mst-1] ) ;
108 }
109
110 } /* end of while stack is non-empty */
111 return ;
112 }
113
114 /*------- sort an array partially recursively, and partially insertion -------*/
115
116 /* quick_sort : sort an array partially recursively, and partially insertion */
117
118 #undef DTYPE
119 #define DTYPE double
120 #include "cs_qsort_small.h"
121
qsort_double(int n,double * a)122 void qsort_double( int n , double *a )
123 {
124 if( n <= 1 ) return ;
125 switch(n){
126 default: break ; /* handled below */
127 case 2: qsort2_double(a) ; return ;
128 case 3: qsort3_double(a) ; return ;
129 case 4: qsort4_double(a) ; return ;
130 case 5: qsort5_double(a) ; return ;
131 case 6: qsort6_double(a) ; return ;
132 case 7: qsort7_double(a) ; return ;
133 case 8: qsort8_double(a) ; return ;
134 case 9: qsort9_double(a) ; return ;
135 case 10: qsort10_double(a); return ;
136 case 11: qsort11_double(a); return ;
137 case 12: qsort12_double(a); return ;
138 case 13: qsort13_double(a); return ;
139 case 14: qsort14_double(a); return ;
140 case 15: qsort15_double(a); return ;
141 case 16: qsort16_double(a); return ;
142 case 17: qsort17_double(a); return ;
143 case 18: qsort18_double(a); return ;
144 case 19: qsort19_double(a); return ;
145 case 20: qsort20_double(a); return ;
146 case 21: qsort21_double(a); return ;
147 case 25: qsort25_double(a); return ;
148 case 27: qsort27_double(a); return ;
149 }
150 qsrec_double( n , a , QS_CUTOFF ) ;
151 isort_double( n , a ) ;
152 return ;
153 }
154