c1-2 (Numerical Recipes in C), страница 2

An immediate consequence is thatthe array elements bb[1], bb[2], bb[3], and bb[4] all exist. In other words therange of bb is bb[1..4]. We will refer to bb as a unit-offset vector. (See AppendixB for some additional discussion of technical details.)It is sometimes convenient to use zero-offset vectors, and sometimes convenientto use unit-offset vectors in algorithms. The choice should be whichever is mostnatural to the problem at hand.

For example, the coefficients of a polynomiala0 + a1 x + a2 x2 + . . . + an xn clearly cry out for the zero-offset a[0..n], whilea vector of N data points xi , i = 1 . . . N calls for a unit-offset x[1..N]. When aroutine in this book has an array as an argument, its header comment always givesthe expected index range. For example,void someroutine(float bb[], int nn)This routine does something with the vector bb[1..nn]....Now, suppose you want someroutine() to do its thing on your own vector,of length 7, say. If your vector, call it aa, is already unit-offset (has the valid rangeaa[1..7]), then you can invoke someroutine(aa,7); in the obvious way.

That isthe recommended procedure, since someroutine() presumably has some logical,or at least aesthetic, reason for wanting a unit-offset vector.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).There is a close, and elegant, correspondence in C between pointers and arrays.The value referenced by an expression like a[j] is defined to be *((a)+(j)),that is, “the contents of the address obtained by incrementing the pointer a byj.” A consequence of this definition is that if a points to a legal data location,the array element a[0] is always defined. Arrays in C are natively “zero-origin”or “zero-offset.” An array declared by the statement float b[4]; has the validreferences b[0], b[1], b[2], and b[3], but not b[4].Right away we need a notation to indicate what is the valid range of an arrayindex.

(The issue comes up about a thousand times in this book!) For the aboveexample, the index range of b will be henceforth denoted b[0..3], a notationborrowed from Pascal. In general, the range of an array declared by floata[M]; is a[0..M − 1], and the same if float is replaced by any other data type.One problem is that many algorithms naturally like to go from 1 to M , notfrom 0 to M − 1.

Sure, you can always convert them, but they then often acquirea baggage of additional arithmetic in array indices that is, at best, distracting. It isbetter to use the power of the C language, in a consistent way, to make the problemdisappear. Consider1.2 Some C Conventions for Scientific Computing19float *vector(long nl, long nh)Allocates a float vector with range [nl..nh].int *ivector(long nl, long nh)Allocates an int vector with range [nl..nh].unsigned char *cvector(long nl, long nh)Allocates an unsigned char vector with range [nl..nh].unsigned long *lvector(long nl, long nh)Allocates an unsigned long vector with range [nl..nh].double *dvector(long nl, long nh)Allocates a double vector with range [nl..nh].A typical use of the above utilities is the declaration float *b; followed byb=vector(1,7);, which makes the range b[1..7] come into existence and allowsb to be passed to any function calling for a unit-offset vector.The file nrutil.c also contains the corresponding deallocation routines,void free_vector(float *v, long nl, long nh)void free_ivector(int *v, long nl, long nh)void free_cvector(unsigned char *v, long nl, long nh)void free_lvector(unsigned long *v, long nl, long nh)void free_dvector(double *v, long nl, long nh)with the typical use being free_vector(b,1,7);.Our recipes use the above utilities extensively for the allocation and deallocationof vector workspace.

We also commend them to you for use in your main programs orother procedures. Note that if you want to allocate vectors of length longer than 64kon an IBM PC-compatible computer, you should replace all occurrences of mallocin nrutil.c by your compiler’s special-purpose memory allocation function. Thisapplies also to matrix allocation, to be discussed next.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).But suppose that your vector of length 7, now call it a, is perversely a native C,zero-offset array (has range a[0..6]). Perhaps this is the case because you disagreewith our aesthetic prejudices, Heaven help you! To use our recipe, do you have tocopy a’s contents element by element into another, unit-offset vector? No! Do youhave to declare a new pointer aaa and set it equal to a-1? No! You simply invokesomeroutine(a-1,7);.

Then a[1], as seen from within our recipe, is actuallya[0] as seen from your program. In other words, you can change conventions “onthe fly” with just a couple of keystrokes.Forgive us for belaboring these points. We want to free you from the zero-offsetthinking that C encourages but (as we see) does not require. A final liberating pointis that the utility file nrutil.c, listed in full in Appendix B, includes functionsfor allocating (using malloc()) arbitrary-offset vectors of arbitrary lengths. Thesynopses of these functions are as follows:20Chapter 1.PreliminariesMatrices and Two-Dimensional Arraysvoid someroutine(a,m,n)float a[m][n];/* ILLEGAL DECLARATION */and emit code to evaluate the variable dimensions m and n (or any variable-dimensionexpression) each time someroutine() is entered.

Alas! the above fragment isforbidden by the C language definition. The implementation of variable dimensionsin C instead requires some additional finesse; however, we will see that one isrewarded for the effort.There is a subtle near-ambiguity in the C syntax for two-dimensional arrayreferences.

Let us elucidate it, and then turn it to our advantage. Consider thearray reference to a (say) float value a[i][j], where i and j are expressionsthat evaluate to type int. A C compiler will emit quite different machine code forthis reference, depending on how the identifier a has been declared. If a has beendeclared as a fixed-size array, e.g., float a[5][9];, then the machine code is: “tothe address a add 9 times i, then add j, return the value thus addressed.” Notice thatthe constant 9 needs to be known in order to effect the calculation, and an integermultiplication is required (see Figure 1.2.1).Suppose, on the other hand, that a has been declared by float **a;.

Thenthe machine code for a[i][j] is: “to the address of a add i, take the value thusaddressed as a new address, add j to it, return the value addressed by this newaddress.” Notice that the underlying size of a[][] does not enter this calculationat all, and that there is no multiplication; an additional indirection replaces it. Wethus have, in general, a faster and more versatile scheme than the previous one.The price that we pay is the storage requirement for one array of pointers (to therows of a[][]), and the slight inconvenience of remembering to initialize thosepointers when we declare an array.Here is our bottom line: We avoid the fixed-size two-dimensional arrays of C asbeing unsuitable data structures for representing matrices in scientific computing.

Weadopt instead the convention “pointer to array of pointers,” with the array elementspointing to the first element in the rows of each matrix. Figure 1.2.1 contrasts therejected and adopted schemes.The following fragment shows how a fixed-size array a of size 13 by 9 isconverted to a “pointer to array of pointers” reference aa:Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The zero- versus unit-offset issue arises here, too.