c21-1 (779629), страница 5
Текст из файла (страница 5)
Besides the declarations for the routines in nrutil.h, itcontains several macros that are used throughout the book. The ideas behind someof these macros were described in §1.2.Memory Allocation: Advanced TopicsTwo issues regarding vector and matrix memory allocation are worth discussing here:first, a “back-door” compatibility between our matrices and those used by other C software;second, the question of whether our pointer arithmetic for unit-offset vectors and matricesviolates the ANSI C standard.1. Back-door compatibility. As explained in §1.2, we always allocate matrices via thescheme “pointer to an array of pointers to rows” (shown in Figure 1.2.1). With this scheme,once storage for a matrix is allocated, the top level pointer (the “name” of the matrix) canbe used to store a matrix of any size whose dimensions are less than or equal to the size ofthe dimensions allocated.
The only time that the allocated size ever again becomes relevantis when the matrix storage is finally deallocated.In this scheme, the rows of a matrix need not be stored in physically contiguous storage;each row has its own pointer addressing it (see Figure 1.2.1). However, the allocation routinesprinted below in fact (and by design) allocate an entire matrix as one contiguous block, viaa single call to malloc(). This “back-door” fact allows you to use matrices created withour routines directly with other software: The address of the first element in a matrix **a(usually &a[1][1] if a has been allocated with a statement like a=matrix(1,m,1,n); butpossibly &a[0][0] if a was created with zero-offsets) is guaranteed to point to the beginningof a contiguous block containing the full physical matrix, stored by rows.
This address can beused as an argument to other software. (Note that this other software will generally also needto know the physical size of the matrix, or at least the number of columns.)2. Unit-offset vectors. In §1.2, we described how a unit-offset vector bb[1..4] could beobtained from a zero-offset array b[0..3] by the pointer arithmetic operation bb=b-1. Sincebb points to one location before b, bb[1] addresses the same element as b[0], and so on.940Sample 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).Appendix B:Utility Routines (nrutil.c)941Here is a listing of nrutil.h:#ifndef _NR_UTILS_H_#define _NR_UTILS_H_static float sqrarg;#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)static double dsqrarg;#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)static double dmaxarg1,dmaxarg2;#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\(dmaxarg1) : (dmaxarg2))static double dminarg1,dminarg2;#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\(dminarg1) : (dminarg2))static float maxarg1,maxarg2;#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\(maxarg1) : (maxarg2))static float minarg1,minarg2;#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\(minarg1) : (minarg2))static long lmaxarg1,lmaxarg2;#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\(lmaxarg1) : (lmaxarg2))static long lminarg1,lminarg2;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).Strictly speaking, this scheme is not blessed by the ANSI C standard. The problemis not the fact that b-1 points to unallocated storage: location b-1 will never be referencedwithout an increment back into the allocated region.
Rather, the problem is that it mighthappen in rare cases (and probably only on a segmented machine) that the expression b-1has no representation at all. If this occurs, then there is no guarantee that the relationb=(b-1)+1 is satisfied.In practice, this is much less of a problem than one might think. We are not aware ofany compiler or machine on which b=(b-n)+n fails to be true for small integer n. Even on asegmented machine, what typically happens is that the compiler stores some (perhaps illegal)representation of b-n, and that b is recovered when n is added back to this representation. Thememory allocation routines in the First Edition of Numerical Recipes in C, in wide use since1988, all have this “problem”, and we have had not even a single report of their failure in thisrespect (notwithstanding the many readers who have told us that theoretically it could fail).We have also communicated to standards bodies the desirability of blessing “b=(b-n)+n” (atleast for some range of n, say n representable as type short) in a future standard, since therewould seem to be no conflict with existing compilers in doing so.Despite the absence of any experimental (as opposed to theoretical) problem, we havetaken some steps in this edition to make our vector and matrix allocation routines more ANSIcompliant.
In the listing that follows, the parameterNR_END is used as a number of extra storagelocations allocated at the beginning of every vector or matrix block, simply for the purposeof making offset pointer references guaranteed-representable.
We set NR_END to a defaultvalue of 1. This has the effect of making all unit-offset allocations (e.g., b=vector(1,7);or a=matrix(1,128,1,128);) be strictly ANSI compliant. With NR_END = 1, the numberof storage locations wasted is fairly negligible. Allocations with offsets other than 1 (e.g.,b=vector(2,10)) are still theoretically non-compliant, but are virtually unknown in ourroutines. If you need to make such allocations, you may wish to consider increasing the valueof NR_END (noting that larger values increase the amount of wasted storage).942Appendix B#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\(lminarg1) : (lminarg2))static int imaxarg1,imaxarg2;#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\(imaxarg1) : (imaxarg2))#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */void nrerror(char error_text[]);float *vector(long nl, long nh);int *ivector(long nl, long nh);unsigned char *cvector(long nl, long nh);unsigned long *lvector(long nl, long nh);double *dvector(long nl, long nh);float **matrix(long nrl, long nrh, long ncl, long nch);double **dmatrix(long nrl, long nrh, long ncl, long nch);int **imatrix(long nrl, long nrh, long ncl, long nch);float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,long newrl, long newcl);float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);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);void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,long ndl, long ndh);#else /* ANSI *//* traditional - K&R */void nrerror();float *vector();Rest of traditional declarations are here on the diskette.#endif /* ANSI */#endif /* _NR_UTILS_H_ */And here is nrutil.c:#include <stdio.h>#include <stddef.h>#include <stdlib.h>#define NR_END 1#define FREE_ARG char*void nrerror(char error_text[])/* Numerical Recipes standard error handler */{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).static int iminarg1,iminarg2;#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\(iminarg1) : (iminarg2))Utility Routines (nrutil.c)943fprintf(stderr,"Numerical Recipes run-time error...\n");fprintf(stderr,"%s\n",error_text);fprintf(stderr,"...now exiting to system...\n");exit(1);}v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));if (!v) nrerror("allocation failure in vector()");return v-nl+NR_END;}int *ivector(long nl, long nh)/* allocate an int vector with subscript range v[nl..nh] */{int *v;v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));if (!v) nrerror("allocation failure in ivector()");return v-nl+NR_END;}unsigned char *cvector(long nl, long nh)/* allocate an unsigned char vector with subscript range v[nl..nh] */{unsigned char *v;v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));if (!v) nrerror("allocation failure in cvector()");return v-nl+NR_END;}unsigned long *lvector(long nl, long nh)/* allocate an unsigned long vector with subscript range v[nl..nh] */{unsigned long *v;v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));if (!v) nrerror("allocation failure in lvector()");return v-nl+NR_END;}double *dvector(long nl, long nh)/* allocate a double vector with subscript range v[nl..nh] */{double *v;v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));if (!v) nrerror("allocation failure in dvector()");return v-nl+NR_END;}float **matrix(long nrl, long nrh, long ncl, long nch)/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */{long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;float **m;/* allocate pointers to rows */m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));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.















