332 Chapter 8.Sorting For"randomly"ordered data,the operations count goes approximately as N1.25,at least for N50,however,Quicksort is generally faster.The program follows: void shell(unsigned long n,float a[]) Sorts an array a]into ascending numerical order by Shell's method(diminishing increment sort).a is replaced on output by its sorted rearrangement.Normally,the argument n should be set to the size of array a,but if n is smaller than this,then only the first n elements of a are sorted.This feature is used in selip. unsigned long i,j,inc; float v; 1nc=1￥ Determine the starting increment. do inc *3; inc++; while (inc a.The process is then repeated on the left and right subarrays independently,and so on. The partitioning process is carried out by selecting some element,say the leftmost,as the partitioning element a.Scan a pointer up the array until you find an element a,and then scan another pointer down from the end of the array until you find an element<a.These two elements are clearly out of place for the final partitioned array,so exchange them.Continue this process until the pointers cross.This is the right place to insert a,and that round of partitioning is done.The
332 Chapter 8. Sorting Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). For “randomly” ordered data, the operations count goes approximately as N 1.25, at least for N 50, however, Quicksort is generally faster. The program follows: void shell(unsigned long n, float a[]) Sorts an array a[] into ascending numerical order by Shell’s method (diminishing increment sort). a is replaced on output by its sorted rearrangement. Normally, the argument n should be set to the size of array a, but if n is smaller than this, then only the first n elements of a are sorted. This feature is used in selip. { unsigned long i,j,inc; float v; inc=1; Determine the starting increment. do { inc *= 3; inc++; } while (inc v) { Inner loop of straight insertion. a[j]=a[j-inc]; j -= inc; if (j 1); } CITED REFERENCES AND FURTHER READING: Knuth, D.E. 1973, Sorting and Searching, vol. 3 of The Art of Computer Programming (Reading, MA: Addison-Wesley), §5.2.1. [1] Sedgewick, R. 1988, Algorithms, 2nd ed. (Reading, MA: Addison-Wesley), Chapter 8. 8.2 Quicksort Quicksort is, on most machines, on average, for large N, the fastest known sorting algorithm. It is a “partition-exchange” sorting method: A “partitioning element” a is selected from the array. Then by pairwise exchanges of elements, the original array is partitioned into two subarrays. At the end of a round of partitioning, the element a is in its final place in the array. All elements in the left subarray are ≤ a, while all elements in the right subarray are ≥ a. The process is then repeated on the left and right subarrays independently, and so on. The partitioning process is carried out by selecting some element, say the leftmost, as the partitioning element a. Scan a pointer up the array until you find an element > a, and then scan another pointer down from the end of the array until you find an element < a. These two elements are clearly out of place for the final partitioned array, so exchange them. Continue this process until the pointers cross. This is the right place to insert a, and that round of partitioning is done. The
8.2 Quicksort 333 question of the best strategy when an element is equal to the partitioning element is subtle;we refer you to Sedgewick [1]for a discussion.(Answer:You should stop and do an exchange. For speed of execution,we do not implement Quicksort using recursion.Thus the algorithm requires an auxiliary array of storage,of length 2 log2 N,which it uses as a push-down stack for keeping track of the pending subarrays.When a subarray has gotten down to some size M,it becomes faster to sort it by straight insertion (88.1),so we will do this.The optimal setting of M is machine dependent,but M =7 is not too far wrong.Some people advocate leaving the short subarrays unsorted until the end,and then doing one giant insertion sort at the end.Since 81 each element moves at most 7 places,this is just as efficient as doing the sorts immediately,and saves on the overhead.However,on modern machines with paged memory,there is increased overhead when dealing with a large array all at once.We have not found any advantage in saving the insertion sorts till the end. As already mentioned,Quicksort's average running time is fast,but its worst case running time can be very slow:For the worst case it is,in fact,an N2 method! And for the most straightforward implementation of Quicksort it turns out that the worst case is achieved for an input array that is already in order!This ordering of the input array might easily occur in practice.One way to avoid this is to use a little random number generator to choose a random element as the partitioning element.Another is to use instead the median of the first.middle.and last elements as Press. of the current subarray. The great speed of Quicksort comes from the simplicity and efficiency of its 9 inner loop.Simply adding one unnecessary test(for example,a test that your pointer has not moved off the end of the array)can almost double the running time!One avoids such unnecessary tests by placing"sentinels"at either end of the subarray being partitioned.The leftmost sentinel is a,the rightmost a.With the 6 "median-of-three"selection of a partitioning element,we can use the two elements that were not the median to be the sentinels for that subarray. Our implementation closely follows [1]: COMPUTING (ISBN 188812920 #include "nrutil.h" #define SWAP(a,b)temp=(a)(a)=(b);(b)=temp; 10.621 #define M 7 #define NSTACK 50 uurggoglrion Here M is the size of subarrays sorted by straight insertion and NSTACK is the required auxiliary Numerical Recipes 43106 storage. void sort(nsigned long n,f1 oat arr门） (outside Sorts an array arr [1..n]into ascending numerical order using the Quicksort algorithm.n is Software. input;arr is replaced on output by its sorted rearrangement. North unsigned long i,ir=n,j,k,1=1,*istack; int jstack=0; float a,tempi istack=lvector(1,NSTACK); for（;;){ Insertion sort when subarray small enough. if(1r-1=1;1--)[ if (arr[i]<a)break; arr[i+1]sarr[i];
8.2 Quicksort 333 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). question of the best strategy when an element is equal to the partitioning element is subtle; we refer you to Sedgewick [1] for a discussion. (Answer: You should stop and do an exchange.) For speed of execution, we do not implement Quicksort using recursion. Thus the algorithm requires an auxiliary array of storage, of length 2 log 2 N, which it uses as a push-down stack for keeping track of the pending subarrays. When a subarray has gotten down to some size M, it becomes faster to sort it by straight insertion (§8.1), so we will do this. The optimal setting of M is machine dependent, but M = 7 is not too far wrong. Some people advocate leaving the short subarrays unsorted until the end, and then doing one giant insertion sort at the end. Since each element moves at most 7 places, this is just as efficient as doing the sorts immediately, and saves on the overhead. However, on modern machines with paged memory, there is increased overhead when dealing with a large array all at once. We have not found any advantage in saving the insertion sorts till the end. As already mentioned, Quicksort’s average running time is fast, but its worst case running time can be very slow: For the worst case it is, in fact, an N 2 method! And for the most straightforward implementation of Quicksort it turns out that the worst case is achieved for an input array that is already in order! This ordering of the input array might easily occur in practice. One way to avoid this is to use a little random number generator to choose a random element as the partitioning element. Another is to use instead the median of the first, middle, and last elements of the current subarray. The great speed of Quicksort comes from the simplicity and efficiency of its inner loop. Simply adding one unnecessary test (for example, a test that your pointer has not moved off the end of the array) can almost double the running time! One avoids such unnecessary tests by placing “sentinels” at either end of the subarray being partitioned. The leftmost sentinel is ≤ a, the rightmost ≥ a. With the “median-of-three” selection of a partitioning element, we can use the two elements that were not the median to be the sentinels for that subarray. Our implementation closely follows [1]: #include "nrutil.h" #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp; #define M 7 #define NSTACK 50 Here M is the size of subarrays sorted by straight insertion and NSTACK is the required auxiliary storage. void sort(unsigned long n, float arr[]) Sorts an array arr[1..n] into ascending numerical order using the Quicksort algorithm. n is input; arr is replaced on output by its sorted rearrangement. { unsigned long i,ir=n,j,k,l=1,*istack; int jstack=0; float a,temp; istack=lvector(1,NSTACK); for (;;) { Insertion sort when subarray small enough. if (ir-l =l;i--) { if (arr[i] <= a) break; arr[i+1]=arr[i];
334 Chapter 8. Sorting arr[i+1]=a; if (jstack ==0)break; ireistack[istack--]; Pop stack and begin a new round of parti- l=istack[jstack--]; tioning. else k=(1+1r)>>1; Choose median of left,center,and right el- SWAP(arr [k],arr[1+1]) ements as partitioning element a.Also if (arr[1]arr[ir]){ rearrange so that a[1]≤a[1+i]≤a[ir]. SWAP(arr[1],arr[ir]) if (arr[1+1]arr[ir]){ SWAP(arr[1+1],arr[ir]) if(arr[1]>arr[1+1])[ granted for 19881992 SWAP(arr[1],arr[1+1]) 1600 1=1+1; Initialize pointers for partitioning j=1r; a=arr[1+1]; Partitioning element. Cambridge for (;;) Beginning of innermost loop. to any from NUMERICAL RECIPES IN do i++;while (arr[i]a); Scan up to find element a. do j--;while (arr[j]a); Scan down to find element a. if (j i)break; (Nort Pointers crossed.Partitioning complete. SWAP(arr[i],arr[j]); Exchange elements. End of innermost loop. America server computer, to make one paper UnN电.t THE arr[1+1]=arr[j]; Insert partitioning element. ART arr[j]=a; istack +2; Push pointers to larger subarray on stack,process smaller subarray immediately. 9 Programs if (jstack NSTACK)nrerror("NSTACK too small in sort."); for their 1f(ir-1+1>=j-1){ istack[istack]=ir; st st istack[jstack-1]=i; 1r=j-1; else to dir istack[jstack]=j-1; istack[jstack-1]=1; rectcustsen 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 1=1; v@cambri 10-621 free_lvector(istack,1,NSTACK); -43108 As usual you can move any other arrays around at the same time as you sort arr.At the risk of being repetitious: (outside North Software. #include "nrutil.h" #define SWAP(a,b)temp=(a)(a)=(b);(b)=temp; Amer #define M 7 #define NSTACK 50 void sort2(unsigned long n,float arr[],float brr[]) Sorts an array arr [1..n]into ascending order using Quicksort,while making the corresponding rearrangement of the array brr[1..n]. unsigned long i,ir=n,j,k,1=1,*istack; int istack=0; float a,b,tempi
334 Chapter 8. Sorting Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). } arr[i+1]=a; } if (jstack == 0) break; ir=istack[jstack--]; Pop stack and begin a new round of partil=istack[jstack--]; tioning. } else { k=(l+ir) >> 1; Choose median of left, center, and right elements as partitioning element a. Also rearrange so that a[l] ≤ a[l+1] ≤ a[ir]. 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; Initialize pointers for partitioning. j=ir; a=arr[l+1]; Partitioning element. for (;;) { Beginning of innermost loop. do i++; while (arr[i] a. do j--; while (arr[j] > a); Scan down to find element 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); } As usual you can move any other arrays around at the same time as you sort arr. At the risk of being repetitious: #include "nrutil.h" #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp; #define M 7 #define NSTACK 50 void sort2(unsigned long n, float arr[], float brr[]) Sorts an array arr[1..n] into ascending order using Quicksort, while making the corresponding rearrangement of the array brr[1..n]. { unsigned long i,ir=n,j,k,l=1,*istack; int jstack=0; float a,b,temp;
8.2 Quicksort 335 istack=lvector(1,NSTACK); for(;;){ Insertion sort when subarray small enough. if(1r-1=1;1-)[ if (arr[i]>1: Choose median of left,center and right el- SWAP(arr [k]arr [1+1]) ements as partitioning element a.Also SWAP(brr [k],brr[1+1]) rearrange so that a[]≤a[1+1]≤a[ir]. (North if (arr[l]arr[ir]){ America to any server computer, SWAP(arr[1],arr[ir]) tusers to make one paper 1988-1992 by Cambridge University Press.Programs THE SWAP(brr[1],brr[ir]) only). ART if (arr[1+1]arr[ir]){ SWAP(arr[1+1],arr[ir]) SWAP(brr[1+1],brr[ir]) copy for their 2 if(arr[1]>arr[1+1])[ strictly prohibited. SWAP(arr[1],arr[1+1]) Copyright (C) SWAP(brr[1],brr[1+1]) 1=1+1; Initialize pointers for partitioning. j=ir; a=arr[1+1]: Partitioning element. b=brr[1+1]; OF SCIENTIFIC COMPUTING(ISBN 0-521 for (;;) Beginning of innermost loop. doi++;while (arr[i]a); Scan up to find element a. v@cambri do j--;while (arr[j]a); Scan down to find element a. if (i=j-1)[ istack[jstack]=ir; istack[istack-1]=i; 1r=j-1; else istack[jstack]=j-1; istack[istack-1]=l; 1=1;
8.2 Quicksort 335 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). istack=lvector(1,NSTACK); for (;;) { Insertion sort when subarray small enough. if (ir-l =l;i--) { if (arr[i] > 1; Choose median of left, center and right elements as partitioning element a. Also rearrange so that a[l] ≤ a[l+1] ≤ a[ir]. SWAP(arr[k],arr[l+1]) SWAP(brr[k],brr[l+1]) if (arr[l] > arr[ir]) { SWAP(arr[l],arr[ir]) SWAP(brr[l],brr[ir]) } if (arr[l+1] > arr[ir]) { SWAP(arr[l+1],arr[ir]) SWAP(brr[l+1],brr[ir]) } if (arr[l] > arr[l+1]) { SWAP(arr[l],arr[l+1]) SWAP(brr[l],brr[l+1]) } i=l+1; Initialize pointers for partitioning. j=ir; a=arr[l+1]; Partitioning element. b=brr[l+1]; for (;;) { Beginning of innermost loop. do i++; while (arr[i] a. do j--; while (arr[j] > a); Scan down to find element NSTACK) nrerror("NSTACK too small in sort2."); 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; }
336 Chapter 8.Sorting 22 You could,in principle,rearrange any number of additional arrays along with brr,but this becomes wasteful as the number of such arrays becomes large.The preferred technique is to make use of an index table,as described in 88.4. CITED REFERENCES AND FURTHER READING: Sedgewick,R.1978,Communications of the ACM,vol.21,pp.847-857.[1] NUMERICAL 8.3 Heapsort While usually not quite as fast as Quicksort,Heapsort is one of our favorite sorting routines.It is a true "in-place"sort,requiring no auxiliary storage.It is an N log2 N process,not only on average,but also for the worst-case order of input data In fact,its worst case is only 20 percent or so worse than its average running time. 3②州 Press. 需 It is beyond our scope to give a complete exposition on the theory of Heapsort. We will mention the general principles,then let you refer to the references [1.21,or analyze the program yourself,if you want to understand the details. A set of N numbers ai,i=1,...,N,is said to form a "heap"if it satisfies the relation a/2≥aj for1≤j/2<j≤N (8.3.1) 61 Here the division in j/2 means "integer divide,"ie.,is an exact integer or else is rounded down to the closest integer.Definition(8.3.1)will make sense if you think of the numbers a;as being arranged in a binary tree,with the top,"boss,"node being a1,the two "underling"nodes being a2 and a3,their four underling nodes being a4 、三 througha7,etc.(See Figure 8.3.1.)In this form,a heap has every"supervisor"greater Numerica 10521 than or equal to its two"supervisees,"down through the levels of the hierarchy. 431 If you have managed to rearrange your array into an order that forms a heap, E Recipes then sorting it is very easy:You pull off the "top of the heap,"which will be the largest element yet unsorted.Then you "promote"to the top of the heap its largest underling.Then you promote its largest underling,and so on.The process is like North what happens(or is supposed to happen)in a large corporation when the chairman of the board retires.You then repeat the whole process by retiring the new chairman of the board.Evidently the whole thing is an N log2 N process,since each retiring chairman leads to log2 N promotions of underlings. Well,how do you arrange the array into a heap in the first place?The answer is again a"sift-up"process like corporate promotion.Imagine that the corporation starts out with N/2 employees on the production line,but with no supervisors.Now a supervisor is hired to supervise two workers.If he is less capable than one of his workers,that one is promoted in his place,and he joins the production line. After supervisors are hired,then supervisors of supervisors are hired,and so on up
336 Chapter 8. Sorting Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). } } } You could, in principle, rearrange any number of additional arrays along with brr, but this becomes wasteful as the number of such arrays becomes large. The preferred technique is to make use of an index table, as described in §8.4. CITED REFERENCES AND FURTHER READING: Sedgewick, R. 1978, Communications of the ACM, vol. 21, pp. 847–857. [1] 8.3 Heapsort While usually not quite as fast as Quicksort, Heapsort is one of our favorite sorting routines. It is a true “in-place” sort, requiring no auxiliary storage. It is an N log2 N process, not only on average, but also for the worst-case order of input data. In fact, its worst case is only 20 percent or so worse than its average running time. It is beyond our scope to give a complete exposition on the theory of Heapsort. We will mention the general principles, then let you refer to the references [1,2], or analyze the program yourself, if you want to understand the details. A set of N numbers ai, i = 1,...,N, is said to form a “heap” if it satisfies the relation aj/2 ≥ aj for 1 ≤ j/2 < j ≤ N (8.3.1) Here the division in j/2 means “integer divide,” i.e., is an exact integer or else is rounded down to the closest integer. Definition (8.3.1) will make sense if you think of the numbers ai as being arranged in a binary tree, with the top, “boss,” node being a1, the two “underling” nodes being a2 and a3, their four underling nodes being a4 through a7, etc. (See Figure 8.3.1.) In this form, a heap has every “supervisor” greater than or equal to its two “supervisees,” down through the levels of the hierarchy. If you have managed to rearrange your array into an order that forms a heap, then sorting it is very easy: You pull off the “top of the heap,” which will be the largest element yet unsorted. Then you “promote” to the top of the heap its largest underling. Then you promote its largest underling, and so on. The process is like what happens (or is supposed to happen) in a large corporation when the chairman of the board retires. You then repeat the whole process by retiring the new chairman of the board. Evidently the whole thing is an N log 2 N process, since each retiring chairman leads to log2 N promotions of underlings. Well, how do you arrange the array into a heap in the first place? The answer is again a “sift-up” process like corporate promotion. Imagine that the corporation starts out with N/2 employees on the production line, but with no supervisors. Now a supervisor is hired to supervise two workers. If he is less capable than one of his workers, that one is promoted in his place, and he joins the production line. After supervisors are hired, then supervisors of supervisors are hired, and so on up
