Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 82
Текст из файла (страница 82)
. . , N . The input quantities n and arr are not changed.{unsigned long i,indxt,ir=n,itemp,j,k,l=1;int jstack=0,*istack;float a;3398.4 Indexing and Rankingoriginalarrayindextableranktablesortedarray5431411821423231556(b)15536144615(a)245831366347227532341326(c)(d)Figure 8.4.1. (a) An unsorted array of six numbers.
(b) Index table, whose entries are pointers tothe elements of (a) in ascending order. (c) Rank table, whose entries are the ranks of the correspondingelements of (a). (d) Sorted array of the elements in (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];340Chapter 8.Sortinga=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);}If you want to sort an array while making the corresponding rearrangement ofseveral or many other arrays, you should first make an index table, then use it torearrange each array in turn.
This requires two arrays of working space: one tohold the index, and another into which an array is temporarily moved, and fromwhich it is redeposited back on itself in the rearranged order. For 3 arrays, theprocedure looks like this:#include "nrutil.h"void sort3(unsigned long n, float ra[], float rb[], float rc[])Sorts an array ra[1..n] into ascending numerical order while making the corresponding rearrangements of the arrays rb[1..n] and rc[1..n].
An index table is constructed via theroutine indexx.{void indexx(unsigned long n, float arr[], unsigned long indx[]);unsigned long j,*iwksp;float *wksp;iwksp=lvector(1,n);wksp=vector(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_vector(wksp,1,n);free_lvector(iwksp,1,n);Make the index table.Save the array ra.Copy it back in rearranged order.Ditto rb.Ditto rc.}The generalization to any other number of arrays is obviously straightforward.A rank table is different from an index table.
A rank table’s jth entry gives the8.5 Selecting the Mth Largest341rank of the jth element of the original array of keys, ranging from 1 (if that elementwas the smallest) to N (if that element was the largest). One can easily constructa rank table from an index table, however:void rank(unsigned long n, unsigned long indx[], unsigned long irank[])Given indx[1..n] as output from the routine indexx, returns an array irank[1..n], thecorresponding table of ranks.{unsigned long j;for (j=1;j<=n;j++) irank[indx[j]]=j;}Figure 8.4.1 summarizes the concepts discussed in this section.8.5 Selecting the Mth LargestSelection is sorting’s austere sister.
(Say that five times quickly!) Where sortingdemands the rearrangement of an entire data array, selection politely asks for a singlereturned value: What is the kth smallest (or, equivalently, the m = N +1−kth largest)element out of N elements? The fastest methods for selection do, unfortunately,rearrange the array for their own computational purposes, typically putting all smallerelements to the left of the kth, all larger elements to the right, and scrambling theorder within each subset. This side effect is at best innocuous, at worst downrightinconvenient. When the array is very long, so that making a scratch copy of it is taxingon memory, or when the computational burden of the selection is a negligible partof a larger calculation, one turns to selection algorithms without side effects, whichleave the original array undisturbed. Such in place selection is slower than the fasterselection methods by a factor of about 10.
We give routines of both types, below.The most common use of selection is in the statistical characterization of a setof data. One often wants to know the median element in an array, or the top andbottom quartile elements. When N is odd, the median is the kth element, withk = (N + 1)/2. When N is even, statistics books define the median as the arithmeticmean of the elements k = N/2 and k = N/2 + 1 (that is, N/2 from the bottomand N/2 from the top). If you accept such pedantry, you must perform two separateselections to find these elements.
For N > 100 we usually define k = N/2 to bethe median element, pedants be damned.The fastest general method for selection, allowing rearrangement, is partitioning, exactly as was done in the Quicksort algorithm (§8.2). Selecting a “random”partition element, one marches through the array, forcing smaller elements to theleft, larger elements to the right.
As in Quicksort, it is important to optimize theinner loop, using “sentinels” (§8.2) to minimize the number of comparisons. Forsorting, one would then proceed to further partition both subsets. For selection,we can ignore one subset and attend only to the one that contains our desired kthelement. Selection by partitioning thus does not need a stack of pending operations,and its operations count scales as N rather than as N log N (see [1]). Comparisonwith sort in §8.2 should make the following routine obvious:342Chapter 8.Sorting#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;float select(unsigned long k, unsigned long n, float arr[])Returns the kth smallest value in the array arr[1..n]. The input array will be rearrangedto have this value in location arr[k], with all smaller elements moved to arr[1..k-1] (inarbitrary order) and all larger elements in arr[k+1..n] (also in arbitrary order).{unsigned long i,ir,j,l,mid;float a,temp;l=1;ir=n;for (;;) {if (ir <= l+1) {Active partition contains 1 or 2 elements.if (ir == l+1 && arr[ir] < arr[l]) {Case of 2 elements.SWAP(arr[l],arr[ir])}return arr[k];} else {mid=(l+ir) >> 1;Choose median of left, center, and right elSWAP(arr[mid],arr[l+1])ements as partitioning element a.
Alsoif (arr[l] > arr[ir]) {rearrange so that arr[l] ≤ arr[l+1],SWAP(arr[l],arr[ir])arr[ir] ≥ arr[l+1].}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);Scan up to find element > a.do j--; while (arr[j] > a);Scan down to find element < a.if (j < i) break;Pointers crossed. Partitioning complete.SWAP(arr[i],arr[j])}End of innermost loop.arr[l+1]=arr[j];Insert partitioning element.arr[j]=a;if (j >= k) ir=j-1;Keep active the partition that contains theif (j <= k) l=i;kth element.}}}In-place, nondestructive, selection is conceptually simple, but it requires a lotof bookkeeping, and it is correspondingly slower.
The general idea is to pick somenumber M of elements at random, to sort them, and then to make a pass throughthe array counting how many elements fall in each of the M + 1 intervals definedby these elements. The kth largest will fall in one such interval — call it the “live”interval. One then does a second round, first picking M random elements in the liveinterval, and then determining which of the new, finer, M + 1 intervals all presentlylive elements fall into. And so on, until the kth element is finally localized within asingle array of size M , at which point direct selection is possible.How shall we pick M ? The number of rounds, logM N = log2 N/ log2 M ,will be smaller if M is larger; but the work to locate each element among M + 1subintervals will be larger, scaling as log2 M for bisection, say.
Each round3438.5 Selecting the Mth Largestrequires looking at all N elements, if only to find those that are still alive, whilethe bisections are dominated by the N that occur in the first round. MinimizingO(N logM N ) + O(N log2 M ) thus yields the resultM ∼2√log2 N(8.5.1)The square root of the logarithm is so slowly varying that secondary considerations ofmachine timing become important. We use M = 64 as a convenient constant value.Two minor additional tricks in the following routine, selip, are (i) augmentingthe set of M random values by an M + 1st, the arithmetic mean, and (ii) choosingthe M random values “on the fly” in a pass through the data, by a method that makeslater values no less likely to be chosen than earlier ones. (The underlying idea is togive element m > M an M/m chance of being brought into the set.
You can proveby induction that this yields the desired result.)#include "nrutil.h"#define M 64#define BIG 1.0e30#define FREEALL free_vector(sel,1,M+2);free_lvector(isel,1,M+2);float selip(unsigned long k, unsigned long n, float arr[])Returns the kth smallest value in the array arr[1..n]. The input array is not altered.{void shell(unsigned long n, float a[]);unsigned long i,j,jl,jm,ju,kk,mm,nlo,nxtmm,*isel;float ahi,alo,sum,*sel;if (k < 1 || k > n || n <= 0) nrerror("bad input to selip");isel=lvector(1,M+2);sel=vector(1,M+2);kk=k;ahi=BIG;alo = -BIG;for (;;) {Main iteration loop, until desired elemm=nlo=0;ment is isolated.sum=0.0;nxtmm=M+1;for (i=1;i<=n;i++) {Make a pass through the whole array.if (arr[i] >= alo && arr[i] <= ahi) {Consider only elements in the current brackets.mm++;if (arr[i] == alo) nlo++;In case of ties for low bracket.Now use statistical procedure for selecting m in-range elements with equalprobability, even without knowing in advance how many there are!if (mm <= M) sel[mm]=arr[i];else if (mm == nxtmm) {nxtmm=mm+mm/M;sel[1 + ((i+mm+kk) % M)]=arr[i]; The % operation provides a some}what random number.sum += arr[i];}}if (kk <= nlo) {Desired element is tied for lower bound;FREEALLreturn it.return alo;}else if (mm <= M) {All in-range elements were kept.
So reshell(mm,sel);turn answer by direct method.ahi = sel[kk];344Chapter 8.SortingFREEALLreturn ahi;}sel[M+1]=sum/mm;Augment selected set by mean value (fixesshell(M+1,sel);degeneracies), and sort it.sel[M+2]=ahi;for (j=1;j<=M+2;j++) isel[j]=0;Zero the count array.for (i=1;i<=n;i++) {Make another pass through the array.if (arr[i] >= alo && arr[i] <= ahi) {For each in-range element..jl=0;ju=M+2;while (ju-jl > 1) {...find its position among the select byjm=(ju+jl)/2;bisection...if (arr[i] >= sel[jm]) jl=jm;else ju=jm;}isel[ju]++;...and increment the counter.}}j=1;Now we can narrow the bounds to justwhile (kk > isel[j]) {one bin, that is, by a factor of orderalo=sel[j];m.kk -= isel[j++];}ahi=sel[j];}}Approximate timings: selip is about 10 times slower than select. Indeed,for N in the range of ∼ 105 , selip is about 1.5 times slower than a full sort withsort, while select is about 6 times faster than sort.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
