1.2 Some C Conventions for Scientific Computing 15 if (julian >IGREG){ Cross-over to Gregorian Calendar produces this correc- ja1pha=(1ong)((doub1e)(julian-1867216)-0.25)/36524.25); tion. ja=julian+1+jalpha-(long)(0.25*jalpha); else if (julian 0){Make day number positive by adding integer number of ja=ju11an+36525*(1-jul1an/36525); Julian centuries,then subtract them off else at the end. ja=julian: jb=ja+1524; jc=(1ong)(6680.0+((doub1e)(jb-2439870)-122.1)/365.25); jd=(1ong)(365*jc+(0.25*jc)); je=(1ong)((jb-jd)/30.6001); *id=jb-jd-(1ong)(30.6001*je); *mm=je-1; 1f(*mm>12)*mm-=12; *1yyy=jc-4715; if(*mm>2)--(*iyyy）; f(*1yyy<=o)-(*1yyy): 1f(ju1ian<0)*1yyy-=100*(1-ju11an/36525): 三今 19881992 11800 NUMERICAL RECIPES I (For additional calendrical algorithms,applicable to various historical calendars,see [8].) CITED REFERENCES AND FURTHER READING: Harbison,S.P,and Steele,G.L.,Jr.1991,C:A Reference Manual,3rd ed.(Englewood Cliffs. NJ:Prentice-Hall). 子 Press. Kernighan,B.W.1978,The Elements of Programming Style (New York:McGraw-Hill).[1] Yourdon,E.1975,Techniques of Program Structure and Design(Englewood Cliffs,NJ:Prentice- Hal0.[2] Jones,R.,and Stewart,I.1987,The Art of C Programming (New York:Springer-Verlag).[3] Hoare,C.A.R.1981,Communications of the ACM,vol.24,pp.75-83. Wirth,N.1983,Programming in Modula-2,3rd ed.(New York:Springer-Verlag).[4] Stroustrup,B.1986,The C++Programming Language(Reading,MA:Addison-Wesley).[5] 6 Borland International,Inc.1989,Turbo Pascal 5.5 Object-Oriented Programming Guide(Scotts Valley,CA:Borland International).[6] Meeus,J.1982,Astronomical Formulae for Calculators,2nd ed.,revised and enlarged(Rich- 令 mond,VA:Willmann-Bell).[7] Hatcher,D.A.1984,Quarterly Journal of the Royal Astronomical Society,vol.25,pp.53-55;see also op.cit1985,vol.26,pp.151-155,and1986,vol.27,Pp.506-507.[8 Numerical Recipes 10621 43106 1.2 Some C Conventions for Scientific (outside Computing North Software. The C language was devised originally for systems programming work,not for scientific computing.Relative to other high-level programming languages,C puts the programmer"very close to the machine"in several respects.It is operator-rich, giving direct access to most capabilities of a machine-language instruction set.It has a large variety of intrinsic data types (short and long,signed and unsigned integers;floating and double-precision reals;pointer types;etc.),and a concise syntax for effecting conversions and indirections.It defines an arithmetic on pointers (addresses)that relates gracefully to array addressing and is highly compatible with the index register structure of many computers
1.2 Some C Conventions for Scientific Computing 15 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). if (julian >= IGREG) { Cross-over to Gregorian Calendar produces this correcjalpha=(long)(((double) (julian-1867216)-0.25)/36524.25); tion. ja=julian+1+jalpha-(long) (0.25*jalpha); } else if (julian 12) *mm -= 12; *iyyy=jc-4715; if (*mm > 2) --(*iyyy); if (*iyyy <= 0) --(*iyyy); if (julian < 0) *iyyy -= 100*(1-julian/36525); } (For additional calendrical algorithms, applicable to various historical calendars, see [8].) CITED REFERENCES AND FURTHER READING: Harbison, S.P., and Steele, G.L., Jr. 1991, C: A Reference Manual, 3rd ed. (Englewood Cliffs, NJ: Prentice-Hall). Kernighan, B.W. 1978, The Elements of Programming Style (New York: McGraw-Hill). [1] Yourdon, E. 1975, Techniques of Program Structure and Design (Englewood Cliffs, NJ: PrenticeHall). [2] Jones, R., and Stewart, I. 1987, The Art of C Programming (New York: Springer-Verlag). [3] Hoare, C.A.R. 1981, Communications of the ACM, vol. 24, pp. 75–83. Wirth, N. 1983, Programming in Modula-2, 3rd ed. (New York: Springer-Verlag). [4] Stroustrup, B. 1986, The C++ Programming Language (Reading, MA: Addison-Wesley). [5] Borland International, Inc. 1989, Turbo Pascal 5.5 Object-Oriented Programming Guide (Scotts Valley, CA: Borland International). [6] Meeus, J. 1982, Astronomical Formulae for Calculators, 2nd ed., revised and enlarged (Richmond, VA: Willmann-Bell). [7] Hatcher, D.A. 1984, Quarterly Journal of the Royal Astronomical Society, vol. 25, pp. 53–55; see also op. cit. 1985, vol. 26, pp. 151–155, and 1986, vol. 27, pp. 506–507. [8] 1.2 Some C Conventions for Scientific Computing The C language was devised originally for systems programming work, not for scientific computing. Relative to other high-level programming languages, C puts the programmer “very close to the machine” in several respects. It is operator-rich, giving direct access to most capabilities of a machine-language instruction set. It has a large variety of intrinsic data types (short and long, signed and unsigned integers; floating and double-precision reals; pointer types; etc.), and a concise syntax for effecting conversions and indirections. It defines an arithmetic on pointers (addresses) that relates gracefully to array addressing and is highly compatible with the index register structure of many computers
16 Chapter 1.Preliminaries Portability has always been another strong point of the C language.C is the underlying language of the UNIX operating system;both the language and the operating system have by now been implemented on literally hundreds of different computers.The language's universality,portability,and flexibility have attracted increasing numbers of scientists and engineers to it.It is commonly used for the real-time control of experimental hardware,often in spite of the fact that the standard UNIX kernel is less than ideal as an operating system for this purpose. The use of C for higher level scientific calculations such as data analysis, modeling,and floating-point numerical work has generally beenslower in developing. 三 In part this is due to the entrenched position of FORTRAN as the mother-tongue of virtually all scientists and engineers born before 1960,and most born after.In part,also,the slowness of C's penetration into scientific computing has been due to deficiencies in the language that computer scientists have been(we think,stubbornly) slow to recognize.Examples are the lack of a good way to raise numbers to small integer powers,and the"implicit conversion of float to double"issue,discussed below.Many,though not all,of these deficiencies are overcome in the ANSI C Standard.Some remaining deficiencies will undoubtedly disappear over time. Yet another inhibition to the mass conversion of scientists to the C cult has been, up to the time of writing,the decided lack of high-quality scientific or numerical 9 libraries.That is the lacuna into which we thrust this edition of Numerical Recipes. We certainly do not claim to be a complete solution to the problem.We do hope to inspire further efforts,and to lay out by example a set of sensible,practical conventions for scientific C programming. The need for programming conventions in C is very great.Far from the problem of overcoming constraints imposed by the language (our repeated experience with OF SCIENTIFIC Pascal),the problem in C is to choose the best and most natural techniques from multiple opportunities-and then to use those techniques completely consistently 6 from program to program.In the rest of this section,we set out some of the issues, and describe the adopted conventions that are used in all of the routines in this book Function Prototypes and Header Files w Numerical Recipes 10621 ANSI C allows functions to be defined with fumnction prototypes,which specify the type of each function parameter.If a function declaration or definition with 43106 a prototype is visible,the compiler can check that a given function call invokes the function with the correct argument types.All the routines printed in this book (outside are in ANSI C prototype form.For the benefit of readers with older "traditional K&R"C compilers,the Numerical Recipes C Diskette includes two complete sets of North Software. programs,one in ANSI,the other in K&R. 言 The easiest way to understand prototypes is by example.A function definition that would be written in traditional c as int g(x,y,z) int x,y; float z; becomes in ANSI C
16 Chapter 1. Preliminaries 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). Portability has always been another strong point of the C language. C is the underlying language of the UNIX operating system; both the language and the operating system have by now been implemented on literally hundreds of different computers. The language’s universality, portability, and flexibility have attracted increasing numbers of scientists and engineers to it. It is commonly used for the real-time control of experimental hardware, often in spite of the fact that the standard UNIX kernel is less than ideal as an operating system for this purpose. The use of C for higher level scientific calculations such as data analysis, modeling, and floating-point numerical work has generally been slower in developing. In part this is due to the entrenched position of FORTRAN as the mother-tongue of virtually all scientists and engineers born before 1960, and most born after. In part, also, the slowness of C’s penetration into scientific computing has been due to deficiencies in the language that computer scientists have been (we think, stubbornly) slow to recognize. Examples are the lack of a good way to raise numbers to small integer powers, and the “implicit conversion of float to double” issue, discussed below. Many, though not all, of these deficiencies are overcome in the ANSI C Standard. Some remaining deficiencies will undoubtedly disappear over time. Yet another inhibition to the mass conversion of scientists to the C cult has been, up to the time of writing, the decided lack of high-quality scientific or numerical libraries. That is the lacuna into which we thrust this edition of Numerical Recipes. We certainly do not claim to be a complete solution to the problem. We do hope to inspire further efforts, and to lay out by example a set of sensible, practical conventions for scientific C programming. The need for programming conventions in C is very great. Far from the problem of overcoming constraints imposed by the language (our repeated experience with Pascal), the problem in C is to choose the best and most natural techniques from multiple opportunities — and then to use those techniques completely consistently from program to program. In the rest of this section, we set out some of the issues, and describe the adopted conventions that are used in all of the routines in this book. Function Prototypes and Header Files ANSI C allows functions to be defined with function prototypes, which specify the type of each function parameter. If a function declaration or definition with a prototype is visible, the compiler can check that a given function call invokes the function with the correct argument types. All the routines printed in this book are in ANSI C prototype form. For the benefit of readers with older “traditional K&R” C compilers, the Numerical Recipes C Diskette includes two complete sets of programs, one in ANSI, the other in K&R. The easiest way to understand prototypes is by example. A function definition that would be written in traditional C as int g(x,y,z) int x,y; float z; becomes in ANSI C
1.2 Some C Conventions for Scientific Computing 17 int g(int x,int y,float z) A function that has no parameters has the parameter type list void. A function declaration (as contrasted to a function definition)is used to "introduce"a function to a routine that is going to call it.The calling routine needs to know the number and type of arguments and the type of the returned value.In a function declaration,you are allowed to omit the parameter names.Thus the declaration for the above function is allowed to be written int g(int,int,float); If a C program consists of multiple source files,the compiler cannot check the consistency of each function call without some additional assistance.The safest 超云 way to proceed is as follows: Every external function should have a single prototype declaration in a header (.h)file. The source file with the definition (body)of the function should also include the header file so that the compiler can check that the prototypes 9 in the declaration and the definition match. Every source file that calls the function should include the appropriate header (.h)file. .Optionally,a routine that calls a function can also include that function's prototype declaration internally.This is often useful when you are 三sn developing a program,since it gives you a visible reminder(checked by the compiler through the common.h file)of a function's argument types. Later,after your program is debugged,you can go back and delete the 的 supernumary internal declarations. For the routines in this book,the header file containing all the prototypes is nr.h, listed in Appendix A.You should put the statement #include nr.h at the top of every source file that contains Numerical Recipes routines.Since,more frequently than not,you will want to include more than one Numerical Recipes routine in a 人念一 10621 single source file,we have not printed this #include statement in front of this book's individual program listings,but you should make sure that it is present in your programs. As backup.and in accordance with the last item on the indented list above. we declare the function prototype of all Numerical Recipes routines that are called 腿 by other Numerical Recipes routines internally to the calling routine.(That also North makes our routines much more readable.)The only exception to this rule is that the small number of utility routines that we use repeatedly (described below)are declared in the additional header file nrutil.h,and the line #include nrutil.h is explicitly printed whenever it is needed. A final important point about the header file nr.h is that,as furnished on the diskette,it contains both ANSI C and traditional K&R-style declarations.The ANSI forms are invoked if any of the following macros are defined:__STDC__ ANSI,or NRANSI.(The purpose of the last name is to give you an invocation that does not conflict with other possible uses of the first two names.If you have an ANSI compiler,it is essential that you invoke it with one or more of these macros
1.2 Some C Conventions for Scientific Computing 17 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). int g(int x, int y, float z) A function that has no parameters has the parameter type list void. A function declaration (as contrasted to a function definition) is used to “introduce” a function to a routine that is going to call it. The calling routine needs to know the number and type of arguments and the type of the returned value. In a function declaration, you are allowed to omit the parameter names. Thus the declaration for the above function is allowed to be written int g(int, int, float); If a C program consists of multiple source files, the compiler cannot check the consistency of each function call without some additional assistance. The safest way to proceed is as follows: • Every external function should have a single prototype declaration in a header (.h) file. • The source file with the definition (body) of the function should also include the header file so that the compiler can check that the prototypes in the declaration and the definition match. • Every source file that calls the function should include the appropriate header (.h) file. • Optionally, a routine that calls a function can also include that function’s prototype declaration internally. This is often useful when you are developing a program, since it gives you a visible reminder (checked by the compiler through the common .h file) of a function’s argument types. Later, after your program is debugged, you can go back and delete the supernumary internal declarations. For the routines in this book, the header file containing all the prototypes is nr.h, listed in Appendix A. You should put the statement #include nr.h at the top of every source file that contains Numerical Recipes routines. Since, more frequently than not, you will want to include more than one Numerical Recipes routine in a single source file, we have not printed this #include statement in front of this book’s individual program listings, but you should make sure that it is present in your programs. As backup, and in accordance with the last item on the indented list above, we declare the function prototype of all Numerical Recipes routines that are called by other Numerical Recipes routines internally to the calling routine. (That also makes our routines much more readable.) The only exception to this rule is that the small number of utility routines that we use repeatedly (described below) are declared in the additional header file nrutil.h, and the line #include nrutil.h is explicitly printed whenever it is needed. A final important point about the header file nr.h is that, as furnished on the diskette, it contains both ANSI C and traditional K&R-style declarations. The ANSI forms are invoked if any of the following macros are defined: __STDC__, ANSI, or NRANSI. (The purpose of the last name is to give you an invocation that does not conflict with other possible uses of the first two names.) If you have an ANSI compiler, it is essential that you invoke it with one or more of these macros
18 Chapter 1.Preliminaries defined.The typical means for doing so is to include a switch like "-DANSI"on the compiler command line. Some further details about the file nr.h are given in Appendix A. Vectors and One-Dimensional Arrays 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 by j."A consequence of this definition is that if a points to a legal data location, the array element a[o]is always defined.Arrays in C are natively "zero-origin" or "zero-offset."An array declared by the statement float b[4];has the valid references b[o].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 array index.(The issue comes up about a thousand times in this book!)For the above example,the index range of b will be henceforth denoted b[o..3],a notation borrowed from Pascal.In general,the range of an array declared by float a[M];is a[0..M-1],and the same if float is replaced by any other data type. 9 One problem is that many algorithms naturally like to go from 1 to M,not from 0 to M-1.Sure,you can always convert them,but they then often acquire a baggage of additional arithmetic in array indices that is,at best,distracting.It is better to use the power of the c language,in a consistent way,to make the problem disappear. Consider S量岂g今 float b[4],*bb; bbsb-1; The pointer bb now points one location before b.An immediate consequence is that the array elements bb[1],bb[2],bb[3],and bb[4]all exist.In other words the range of bb is bb[1..4].We will refer to bb as a unit-offset vector.(See Appendix B for some additional discussion of technical details.) It is sometimes convenient to use zero-offset vectors,and sometimes convenient 10.621 to use unit-offset vectors in algorithms.The choice should be whichever is most Numerica natural to the problem at hand.For example,the coefficients of a polynomial 4310 ao +ax+a2x2+...+anz"clearly cry out for the zero-offset a[o..n],while Recipes a vector of N data pointsi,i=1...N calls for a unit-offset x[1..N].When a routine in this book has an array as an argument,its header comment always gives the expected index range.For example, North Software. 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 range aa [1..7]).then you can invoke someroutine (aa,7);in the obvious way.That is the recommended procedure,since someroutine (presumably has some logical, or at least aesthetic,reason for wanting a unit-offset vector
18 Chapter 1. Preliminaries 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). defined. The typical means for doing so is to include a switch like “-DANSI” on the compiler command line. Some further details about the file nr.h are given in Appendix A. Vectors and One-Dimensional Arrays 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 by j.” 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 valid references 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 array index. (The issue comes up about a thousand times in this book!) For the above example, the index range of b will be henceforth denoted b[0..3], a notation borrowed from Pascal. In general, the range of an array declared by float a[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, not from 0 to M − 1. Sure, you can always convert them, but they then often acquire a baggage of additional arithmetic in array indices that is, at best, distracting. It is better to use the power of the C language, in a consistent way, to make the problem disappear. Consider float b[4],*bb; bb=b-1; The pointer bb now points one location before b. An immediate consequence is that the array elements bb[1], bb[2], bb[3], and bb[4] all exist. In other words the range of bb is bb[1..4]. We will refer to bb as a unit-offset vector. (See Appendix B for some additional discussion of technical details.) It is sometimes convenient to use zero-offset vectors, and sometimes convenient to use unit-offset vectors in algorithms. The choice should be whichever is most natural to the problem at hand. For example, the coefficients of a polynomial a0 + a1x + a2x2 + ... + anxn clearly cry out for the zero-offset a[0..n], while a vector of N data points xi, i = 1 ...N calls for a unit-offset x[1..N]. When a routine in this book has an array as an argument, its header comment always gives the 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 range aa[1..7]), then you can invoke someroutine(aa,7); in the obvious way. That is the 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 Computing 19 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 disagree with our aesthetic prejudices,Heaven help you!To use our recipe,do you have to copy a's contents element by element into another,unit-offset vector?No!Do you have to declare a new pointer aaa and set it equal to a-1?No!You simply invoke someroutine(a-1,7);.Then a[1],as seen from within our recipe,is actually a[o]as seen from your program.In other words,you can change conventions"on the fly"with just a couple of keystrokes. Forgive us for belaboring these points.We want to free you from the zero-offset thinking that C encourages but(as we see)does not require.A final liberating point is that the utility file nrutil.c,listed in full in Appendix B,includes functions for allocating (using malloc())arbitrary-offset vectors of arbitrary lengths.The nted for synopses of these functions are as follows: float *vector(long nl,long nh) Allocates a float vector with range [nl..nh] to any 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]. (North America server computer, to make one paper UnN电.t THE unsigned long *lvector(long nl,long nh) ART Allocates an unsigned long vector with range [nl..nh]. Programs double *dvector(long nl,long nh) Allocates a double vector with range [nl..nh]. strictly proh A typical use of the above utilities is the declaration float *b;followed by b=vector(1,7);,which makes the range b[1..7]come into existence and allows b to be passed to any function calling for a unit-offset vector. The file nrutil.c also contains the corresponding deallocation routines, 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN void free_vector(float *v,long nl,long nh) void free_ivector(int *v,long nl,long nh) 10-621 void free_cvector(unsigned char *v,long nl,long nh) 43108 void free_lvector(unsigned long *v,long nl,long nh) void free_dvector(double *v,long nl,long nh) (outside North Software. with the typical use being free_vector(b,1,7);. Our recipes use the above utilities extensively for the allocation and deallocation Ame of vector workspace.We also commend them to you for use in your main programs or other procedures.Note that if you want to allocate vectors of length longer than 64k on an IBM PC-compatible computer,you should replace all occurrences of malloc in nrutil.c by your compiler's special-purpose memory allocation function.This applies also to matrix allocation,to be discussed next
1.2 Some C Conventions for Scientific Computing 19 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). 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 disagree with our aesthetic prejudices, Heaven help you! To use our recipe, do you have to copy a’s contents element by element into another, unit-offset vector? No! Do you have to declare a new pointer aaa and set it equal to a-1? No! You simply invoke someroutine(a-1,7);. Then a[1], as seen from within our recipe, is actually a[0] as seen from your program. In other words, you can change conventions “on the fly” with just a couple of keystrokes. Forgive us for belaboring these points. We want to free you from the zero-offset thinking that C encourages but (as we see) does not require. A final liberating point is that the utility file nrutil.c, listed in full in Appendix B, includes functions for allocating (using malloc()) arbitrary-offset vectors of arbitrary lengths. The synopses 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 by b=vector(1,7);, which makes the range b[1..7] come into existence and allows b 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 deallocation of vector workspace. We also commend them to you for use in your main programs or other procedures. Note that if you want to allocate vectors of length longer than 64k on an IBM PC-compatible computer, you should replace all occurrences of malloc in nrutil.c by your compiler’s special-purpose memory allocation function. This applies also to matrix allocation, to be discussed next
20 Chapter 1.Preliminaries Matrices and Two-Dimensional Arrays The zero-versus unit-offset issue arises here,too.Let us,however,defer it for a moment in favor of an even more fundamental matter.that of variable dimension arrays(FORTRAN terminology)or conformant arrays(Pascal terminology).These are arrays that need to be passed to a function along with real-time information about their two-dimensional size.The systems programmer rarely deals with two- dimensional arrays,and almost never deals with two-dimensional arrays whose size is variable and known only at run time.Such arrays are,however,the bread and butter of scientific computing.Imagine trying to live with a matrix inversion routine that could work with only one size of matrix! There is no technical reason that a C compiler could not allow a syntax like void someroutine(a,m,n) float a[m][n]; /ILLEGAL DECLARATION * y and emit code to evaluate the variable dimensions m and n(or any variable-dimension 令 expression)each time someroutine (is entered.Alas!the above fragment is forbidden by the C language definition.The implementation of variable dimensions in C instead requires some additional finesse;however,we will see that one is Press. rewarded for the effort. There is a subtle near-ambiguity in the C syntax for two-dimensional array 9 references.Let us elucidate it,and then turn it to our advantage.Consider the array reference to a(say)float value a[i][j],where i and j are expressions that evaluate to type int.A C compiler will emit quite different machine code for OF SCIENTIFIC( this reference,depending on how the identifier a has been declared.If a has been 6 declared as a fixed-size array,e.g.,float a[5][9];,then the machine code is:"to the address a add 9 times i.then add i.return the value thus addressed."Notice that the constant 9 needs to be known in order to effect the calculation,and an integer multiplication is required (see Figure 1.2.1). Suppose,on the other hand,that a has been declared by float **a;.Then the machine code for a[i]j]is:"to the address of a add i,take the value thus Numerica 10621 addressed as a new address,add i to it.return the value addressed by this new address."Notice that the underlying size of a]]does not enter this calculation at all,and that there is no multiplication;an additional indirection replaces it.We thus 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 the rows of a[][]),and the slight inconvenience of remembering to initialize those pointers when we declare an array. Here is our bottom line:We avoid the fixed-size two-dimensional arrays of C as being unsuitable data structures for representing matrices in scientific computing.We adopt instead the convention"pointer to array of pointers,"with the array elements pointing to the first element in the rows of each matrix.Figure 1.2.1 contrasts the rejected and adopted schemes. The following fragment shows how a fixed-size array a of size 13 by 9 is converted to a"pointer to array of pointers"reference aa:
20 Chapter 1. Preliminaries 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). Matrices and Two-Dimensional Arrays The zero- versus unit-offset issue arises here, too. Let us, however, defer it for a moment in favor of an even more fundamental matter, that of variable dimension arrays (FORTRAN terminology) or conformant arrays (Pascal terminology). These are arrays that need to be passed to a function along with real-time information about their two-dimensional size. The systems programmer rarely deals with twodimensional arrays, and almost never deals with two-dimensional arrays whose size is variable and known only at run time. Such arrays are, however, the bread and butter of scientific computing. Imagine trying to live with a matrix inversion routine that could work with only one size of matrix! There is no technical reason that a C compiler could not allow a syntax like void someroutine(a,m,n) float a[m][n]; /* ILLEGAL DECLARATION */ and emit code to evaluate the variable dimensions m and n (or any variable-dimension expression) each time someroutine() is entered. Alas! the above fragment is forbidden by the C language definition. The implementation of variable dimensions in C instead requires some additional finesse; however, we will see that one is rewarded for the effort. There is a subtle near-ambiguity in the C syntax for two-dimensional array references. Let us elucidate it, and then turn it to our advantage. Consider the array reference to a (say) float value a[i][j], where i and j are expressions that evaluate to type int. A C compiler will emit quite different machine code for this reference, depending on how the identifier a has been declared. If a has been declared as a fixed-size array, e.g., float a[5][9];, then the machine code is: “to the address a add 9 times i, then add j, return the value thus addressed.” Notice that the constant 9 needs to be known in order to effect the calculation, and an integer multiplication is required (see Figure 1.2.1). Suppose, on the other hand, that a has been declared by float **a;. Then the machine code for a[i][j] is: “to the address of a add i, take the value thus addressed as a new address, add j to it, return the value addressed by this new address.” Notice that the underlying size of a[][] does not enter this calculation at all, and that there is no multiplication; an additional indirection replaces it. We thus 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 the rows of a[][]), and the slight inconvenience of remembering to initialize those pointers when we declare an array. Here is our bottom line: We avoid the fixed-size two-dimensional arrays of C as being unsuitable data structures for representing matrices in scientific computing. We adopt instead the convention “pointer to array of pointers,” with the array elements pointing to the first element in the rows of each matrix. Figure 1.2.1 contrasts the rejected and adopted schemes. The following fragment shows how a fixed-size array a of size 13 by 9 is converted to a “pointer to array of pointers” reference aa:
1.2 Some C Conventions for Scientific Computing 21 [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] (a) http://www.nr. 83g *m[0] [0][0 [0][1] [0][2] 0][3] [0][4] granted for 100 (including this one) 19881992 -872 *m[1] [1][0] [1][1] [1][2] 1][3] [1][4] to any from NUMERICAL RECIPES IN by Cambridge (Nort *m[2] 2][0 2][1 2][2 2][3] 2][4] (b) America server computer, users to make one paper UnN电.t THE Figure 1.2.1.Two storage schemes for a matrix m.Dotted lines denote address reference, while solid ART lines connect sequential memory locations.(a)Pointer to a fixed size two-dimensional array.(b)Pointer to an array of pointers to rows;this is the scheme adopted in this book 9 Programs 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][o] to dir The identifier aa is now a matrix with index range aa [0..12][0..8].You can use OF SCIENTIFIC COMPUTING (ISBN or modify its elements ad lib,and more importantly you can pass it as an argument 1988.1820 to any function by its name aa.That function,which declares the corresponding dummy argument as float **aa;,can address its elements as aa [i][i]without v@cam knowing its physical size. .Further reproduction, 10-621 You may rightly not wish to clutter your programs with code like the above fragment.Also,there is still the outstanding problem of how to treat unit-offset Numerical Recipes 431085 indices,so that(for example)the above matrix aa could be addressed with the range a[1..13][1..9].Both of these problems are solved by additional utility routines (outside in nrutil.c(Appendix B)which allocate and deallocate matrices of arbitrary range. North Software. The synopses are Ame float **matrix(long nrl,long nrh,long ncl,long nch) Allocates a float matrix with range [nrl..nrh][ncl..nch]. double **dmatrix(long nrl,long nrh,long ncl,long nch) Allocates a double matrix with range [nrl..nrh][ncl..nch]. int **imatrix(long nrl,long nrh,long ncl,long nch) Allocates an int matrix with range [nrl..nrh][ncl..nch]. void free_matrix(float **m,long nrl,long nrh,long ncl,long nch) Frees a matrix allocated with matrix
1.2 Some C Conventions for Scientific Computing 21 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). [0][0] [0][1] [0][2] [0][3] [0][4] [1][0] [1][1] [1][2] [1][3] [1][4] [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] [2][0] [2][1] [2][2] [2][3] [2][4] *m[0] *m[1] *m[2] **m **m (a) (b) Figure 1.2.1. Two storage schemes for a matrix m. Dotted lines denote address reference, while solid lines connect sequential memory locations. (a) Pointer to a fixed size two-dimensional array. (b) Pointer to 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 use or modify its elements ad lib, and more importantly you can pass it as an argument to any function by its name aa. That function, which declares the corresponding dummy argument as float **aa;, can address its elements as aa[i][j] without knowing its physical size. You may rightly not wish to clutter your programs with code like the above fragment. Also, there is still the outstanding problem of how to treat unit-offset indices, so that (for example) the above matrix aa could be addressed with the range a[1..13][1..9]. Both of these problems are solved by additional utility routines in nrutil.c (Appendix B) which allocate and deallocate matrices of arbitrary range. The synopses are float **matrix(long nrl, long nrh, long ncl, long nch) Allocates a float matrix with range [nrl..nrh][ncl..nch]. double **dmatrix(long nrl, long nrh, long ncl, long nch) Allocates a double matrix with range [nrl..nrh][ncl..nch]. int **imatrix(long nrl, long nrh, long ncl, long nch) Allocates an int matrix with range [nrl..nrh][ncl..nch]. void free_matrix(float **m, long nrl, long nrh, long ncl, long nch) Frees a matrix allocated with matrix
22 Chapter 1.Preliminaries void free_dmatrix(double **m,long nrl,long nrh,long ncl,long nch) Frees a matrix allocated with dmatrix. void free_imatrix(int **m,long nrl,long nrh,long ncl,long nch) Frees a matrix allocated with imatrix. A typical use is float **a; http://www.nr. Permission is readable files a=matrix(1,13,1,9); a[3][5]=... 83g ..+a[2][9]/3.0. someroutine(a,...); 11-800 (including this one) 19881992 free_matrix(a,1,13,1,9); All matrices in Numerical Recipes are handled with the above paradigm,and we commend it to you. Some further utilities for handling matrices are also included in nrutil.c. (North The first is a function submatrix()that sets up a new pointer reference to an America by Cambridge University Press. users to make one paper from NUMERICAL RECIPES IN computer, THE already-existing matrix (or sub-block thereof),along with new offsets if desired. Its synopsis is ART 9 Programs float **submatrix(float **a,long oldrl,long oldrh,long oldcl, long oldch,long newrl,long nevcl) Point a submatrix [newrl..newrl+(oldrh-oldrl)][newcl..newcl+(oldch-oldcl)]to the existing matrix range a[oldr1..oldrh][oldcl..oldch]. 兰三会袋 to dir Here oldrl and oldrh are respectively the lower and upper row indices of the original matrix that are to be represented by the new matrix,oldcl and oldch are OF SCIENTIFIC COMPUTING (ISBN the corresponding column indices,and newrl and newcl are the lower row and 19841820 column indices for the new matrix.(We don't need upper row and column indices, since they are implied by the quantities already given.) 10-521 Two sample uses might be,first,to select as a 2 x 2 submatrix b[1..2] [1..2]some interior range of an existing matrix,say a[4..5][2..3], Numerical Recipes 43108 f10at**a,**b; (outside a=matrix(1,13,1,9); North Software. b=submatrix(a,4,5,2,3,1,1); and second,to map an existing matrix a[1..13][1..9]into a new matrix b[0.12][0..8], float **a,**b; a=matrix(1,13,1,9): b=submatrix(a,1,13,1,9,0,0）;
22 Chapter 1. Preliminaries 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). void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch) Frees a matrix allocated with dmatrix. void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch) Frees a matrix allocated with imatrix. A typical use is float **a; a=matrix(1,13,1,9); ... a[3][5]=... ...+a[2][9]/3.0... someroutine(a,...); ... free_matrix(a,1,13,1,9); All matrices in Numerical Recipes are handled with the above paradigm, and we commend it to you. Some further utilities for handling matrices are also included in nrutil.c. The first is a function submatrix() that sets up a new pointer reference to an already-existing matrix (or sub-block thereof), along with new offsets if desired. Its synopsis is float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, long newrl, long newcl) Point a submatrix [newrl..newrl+(oldrh-oldrl)][newcl..newcl+(oldch-oldcl)] to the existing matrix range a[oldrl..oldrh][oldcl..oldch]. Here oldrl and oldrh are respectively the lower and upper row indices of the original matrix that are to be represented by the new matrix, oldcl and oldch are the corresponding column indices, and newrl and newcl are the lower row and column indices for the new matrix. (We don’t need upper row and column indices, since they are implied by the quantities already given.) Two sample uses might be, first, to select as a 2 × 2 submatrix b[1..2] [1..2] some interior range of an existing matrix, say a[4..5][2..3], float **a,**b; a=matrix(1,13,1,9); ... b=submatrix(a,4,5,2,3,1,1); and second, to map an existing matrix a[1..13][1..9] into a new matrix b[0..12][0..8], float **a,**b; a=matrix(1,13,1,9); ... b=submatrix(a,1,13,1,9,0,0);
1.2 Some C Conventions for Scientific Computing 23 Incidentally,you can use submatrix()for matrices of any type whose sizeof() is the same as sizeof(float)(often true for int,e.g.);just cast the first argument to type float *and cast the result to the desired type,e.g.,int ** The function void free_submatrix(float **b,long nrl,long nrh,long ncl,long nch) frees the array of row-pointers allocated by submatrix().Note that it does not free the memory allocated to the data in the submatrix,since that space still lies within the memory allocation of some original matrix. 81 Finally,if you have a standard C matrix declared as a[nrow][ncol],and you want to convert it into a matrix declared in our pointer-to-row-of-pointers manner, the following function does the trick: float **convert_matrix(float *a,long nrl,long nrh,long ncl,long nch) Allocate a float matrix m[nrl..nrh][ncl..nch]that points to the matrix declared in the standard C manner as a[nrow][ncol],where nrow=nrh-nr1+1 and ncol-nch-ncl+1.The routine should be called with the address &a[o][o]as the first argument. 令 (You can use this function when you want to make use of C's initializer syntax 造 Press. to set values for a matrix,but then be able to pass the matrix to programs in this book.)The function void free_convert_matrix(float **b,long nrl,long nrh,long ncl,long nch) Free a matrix allocated by convert_matrix() SCIENTIFIC 6 frees the allocation,without affecting the original matrix a. The only examples of allocating a three-dimensional array as a pointer-to- pointer-to-pointer structure in this book are found in the routines rlft3 in 812.5 and sfroid in $17.4.The necessary allocation and deallocation functions are Numerical Recipes 10-621 float ***f3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh) Allocate a float 3-dimensional array with subscript range [nrl..nrh][ncl..nch][ndl..ndh]. uction 43106 void free_f3tensor(float ***t,long nrl,long nrh,long ncl,long nch, long ndl,long ndh) (outside Free a float 3-dimensional array allocated by f3tensor(). Software. Complex Arithmetic C does not have complex data types,or predefined arithmetic operations on complex numbers.That omission is easily remedied with the set of functions in the file complex.c which is printed in full in Appendix C at the back of the book. A synopsis is as follows:
1.2 Some C Conventions for Scientific Computing 23 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). Incidentally, you can use submatrix() for matrices of any type whose sizeof() is the same as sizeof(float) (often true for int, e.g.); just cast the first argument to type float ** and cast the result to the desired type, e.g., int **. The function void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch) frees the array of row-pointers allocated by submatrix(). Note that it does not free the memory allocated to the data in the submatrix, since that space still lies within the memory allocation of some original matrix. Finally, if you have a standard C matrix declared as a[nrow][ncol], and you want to convert it into a matrix declared in our pointer-to-row-of-pointers manner, the following function does the trick: float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch) Allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 and ncol=nch-ncl+1. The routine should be called with the address &a[0][0] as the first argument. (You can use this function when you want to make use of C’s initializer syntax to set values for a matrix, but then be able to pass the matrix to programs in this book.) The function void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch) Free a matrix allocated by convert_matrix(). frees the allocation, without affecting the original matrix a. The only examples of allocating a three-dimensional array as a pointer-topointer-to-pointer structure in this book are found in the routines rlft3 in §12.5 and sfroid in §17.4. The necessary allocation and deallocation functions are float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh) Allocate a float 3-dimensional array with subscript range [nrl..nrh][ncl..nch][ndl..ndh]. void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh) Free a float 3-dimensional array allocated by f3tensor(). Complex Arithmetic C does not have complex data types, or predefined arithmetic operations on complex numbers. That omission is easily remedied with the set of functions in the file complex.c which is printed in full in Appendix C at the back of the book. A synopsis is as follows:
24 Chapter 1.Preliminaries typedef struct FCOMPLEX {float r,i;}fcomplex; fcomplex Cadd(fcomplex a,fcomplex b) Returns the complex sum of two complex numbers fcomplex Csub(fcomplex a,fcomplex b) Returns the complex difference of two complex numbers fcomplex Cmul(fcomplex a,fcomplex b) Returns the complex product of two complex numbers. fcomplex Cdiv(fcomplex a,fcomplex b) read able files Returns the complex quotient of two complex numbers. fcomplex Csqrt(fcomplex z) 83g Returns the complex square root of a complex number. granted for fcomplex Conjg(fcomplex z) 11-800-872 (including this one) Returns the complex conjugate of a complex number float Cabs(fcomplex z) Returns the absolute value (modulus)of a complex number -7423 (North America to any server computer, by Cambridge University Press. users to make one paper from NUMERICAL RECIPES IN C: 198891992 fcomplex Complex(float re,float im) Returns a complex number with specified real and imaginary parts. THE fcomplex RCmul(float x,fcomplex a) Returns the complex product of a real number and a complex number. ART The implementation of several of these complex operations in floating-point arithmetic is not entirely trivial;see $5.4. ictly proh Programs Only about halfa dozen routines in this book make explicit use of these complex arithmetic functions.The resulting code is not as readable as one would like.because the familiar operations +-*are replaced by function calls.The C++extension to to dir the C language allows operators to be redefined.That would allow more readable code.However.in this book we are committed to standard C. OF SCIENTIFIC COMPUTING (ISBN We should mention that the above functions assume the ability to pass,return, 1988-1992 and assign structures like FCOMPLEX(or types such as fcomplex that are defined to be structures)by value.All recent C compilers have this ability,but it is not in Fuurggoglrion Numerical Recipes 10-621 the original K&R C definition.If you are missing it,you will have to rewrite the functions in complex.c,making them pass and return pointers to variables of type 43108 fcomplex instead of the variables themselves.Likewise,you will need to modify the recipes that use the functions. Several other routines (e.g.,the Fourier transforms four1 and fourn)do (outside complex arithmetic"by hand,"that is,they carry around real and imaginary parts as North Software. float variables.This results in more efficient code than would be obtained by using the functions in complex.c.But the code is even less readable.There is simply no ideal solution to the complex arithmetic problem in C. Implicit Conversion of Float to Double In traditional(K&R)C,float variables are automatically converted to double before any operation is attempted,including both arithmetic operations and passing as arguments to functions.All arithmetic is then done in double precision.If a float variable receives the result of such an arithmetic operation,the high precision
24 Chapter 1. Preliminaries 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). typedef struct FCOMPLEX {float r,i;} fcomplex; fcomplex Cadd(fcomplex a, fcomplex b) Returns the complex sum of two complex numbers. fcomplex Csub(fcomplex a, fcomplex b) Returns the complex difference of two complex numbers. fcomplex Cmul(fcomplex a, fcomplex b) Returns the complex product of two complex numbers. fcomplex Cdiv(fcomplex a, fcomplex b) Returns the complex quotient of two complex numbers. fcomplex Csqrt(fcomplex z) Returns the complex square root of a complex number. fcomplex Conjg(fcomplex z) Returns the complex conjugate of a complex number. float Cabs(fcomplex z) Returns the absolute value (modulus) of a complex number. fcomplex Complex(float re, float im) Returns a complex number with specified real and imaginary parts. fcomplex RCmul(float x, fcomplex a) Returns the complex product of a real number and a complex number. The implementation of several of these complex operations in floating-point arithmetic is not entirely trivial; see §5.4. Only about half a dozen routines in this book make explicit use of these complex arithmetic functions. The resulting code is not as readable as one would like, because the familiar operations +-*/ are replaced by function calls. The C++ extension to the C language allows operators to be redefined. That would allow more readable code. However, in this book we are committed to standard C. We should mention that the above functions assume the ability to pass, return, and assign structures like FCOMPLEX (or types such as fcomplex that are defined to be structures) by value. All recent C compilers have this ability, but it is not in the original K&R C definition. If you are missing it, you will have to rewrite the functions in complex.c, making them pass and return pointers to variables of type fcomplex instead of the variables themselves. Likewise, you will need to modify the recipes that use the functions. Several other routines (e.g., the Fourier transforms four1 and fourn) do complex arithmetic “by hand,” that is, they carry around real and imaginary parts as float variables. This results in more efficient code than would be obtained by using the functions in complex.c. But the code is even less readable. There is simply no ideal solution to the complex arithmetic problem in C. Implicit Conversion of Float to Double In traditional (K&R) C, float variables are automatically converted to double before any operation is attempted, including both arithmetic operations and passing as arguments to functions. All arithmetic is then done in double precision. If a float variable receives the result of such an arithmetic operation, the high precision