Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C, страница 10
Описание файла
PDF-файл из архива "Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C", который расположен в категории "". Всё это находится в предмете "численные методы" из 6 семестр, которые можно найти в файловом архиве . Не смотря на прямую связь этого архива с , его также можно найти и в других разделах. Архив можно найти в разделе "книги и методические указания", в предмете "численные методы и алгоритмы" в общих файлах.
Просмотр PDF-файла онлайн
Текст 10 страницы из PDF
(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. Considerfloat b[4],*bb;bb=b-1;The pointer bb now points one location before b. 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.1.2 Some C Conventions for Scientific Computing19But 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:float *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.20Chapter 1.PreliminariesMatrices and Two-Dimensional ArraysThe zero- versus unit-offset issue arises here, too. Let us, however, defer it fora moment in favor of an even more fundamental matter, that of variable dimensionarrays (FORTRAN terminology) or conformant arrays (Pascal terminology). Theseare arrays that need to be passed to a function along with real-time informationabout their two-dimensional size. The systems programmer rarely deals with twodimensional arrays, and almost never deals with two-dimensional arrays whose sizeis variable and known only at run time.
Such arrays are, however, the bread andbutter of scientific computing. Imagine trying to live with a matrix inversion routinethat could work with only one size of matrix!There is no technical reason that a C compiler could not allow a syntax likevoid 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:211.2 Some C Conventions for Scientific Computing**m[0][0][0][1][0][2][0][3][0][4][1][0][1][1][1][2][1][3][1][4][2][0][2][1][2][2][2][3][2][4]*m[0][0][0][0][1][0][2][0][3][0][4]*m[1][1][0][1][1][1][2][1][3][1][4]*m[2][2][0][2][1][2][2][2][3][2][4](a)**m(b)Figure 1.2.1.
Two storage schemes for a matrix m. Dotted lines denote address reference, while solidlines connect sequential memory locations. (a) Pointer to a fixed size two-dimensional array. (b) Pointerto an array of pointers to rows; this is the scheme adopted in this book.float a[13][9],**aa;int i;aa=(float **) malloc((unsigned) 13*sizeof(float*));for(i=0;i<=12;i++) aa[i]=a[i];a[i] is a pointer to a[i][0]The identifier aa is now a matrix with index range aa[0..12][0..8]. You can useor modify its elements ad lib, and more importantly you can pass it as an argumentto any function by its name aa.