c1-0 (Numerical Recipes in C)

Preliminaries1.0 IntroductionThis book, like its predecessor edition, is supposed to teach you methods ofnumerical computing that are practical, efficient, and (insofar as possible) elegant.We presume throughout this book that you, the reader, have particular tasks that youwant to get done.

We view our job as educating you on how to proceed. Occasionallywe may try to reroute you briefly onto a particularly beautiful side road; but by andlarge, we will guide you along main highways that lead to practical destinations.Throughout this book, you will find us fearlessly editorializing, telling youwhat you should and shouldn’t do. This prescriptive tone results from a consciousdecision on our part, and we hope that you will not find it irritating. We do notclaim that our advice is infallible! Rather, we are reacting against a tendency, inthe textbook literature of computation, to discuss every possible method that hasever been invented, without ever offering a practical judgment on relative merit.

Wedo, therefore, offer you our practical judgments whenever we can. As you gainexperience, you will form your own opinion of how reliable our advice is.We presume that you are able to read computer programs in C, that beingthe language of this version of Numerical Recipes (Second Edition). The bookNumerical Recipes in FORTRAN (Second Edition) is separately available, if youprefer to program in that language.

Earlier editions of Numerical Recipes in Pascaland Numerical Recipes Routines and Examples in BASIC are also available; whilenot containing the additional material of the Second Edition versions in C andFORTRAN, these versions are perfectly serviceable if Pascal or BASIC is yourlanguage of choice.When we include programs in the text, they look like this:#include <math.h>#define RAD (3.14159265/180.0)void flmoon(int n, int nph, long *jd, float *frac)Our programs begin with an introductory comment summarizing their purpose and explainingtheir calling sequence.

This routine calculates the phases of the moon. Given an integer n anda code nph for the phase desired (nph = 0 for new moon, 1 for first quarter, 2 for full, 3 for lastquarter), the routine returns the Julian Day Number jd, and the fractional part of a day fracto be added to it, of the nth such phase since January, 1900. Greenwich Mean Time is assumed.{void nrerror(char error_text[]);int i;float am,as,c,t,t2,xtra;This is how we comment an individualline.c=n+nph/4.0;1Sample 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).Chapter 1.2Chapter 1.Preliminaries}If the syntax of the function definition above looks strange to you, then you areprobably used to the older Kernighan and Ritchie (“K&R”) syntax, rather than that ofthe newer ANSI C. In this edition, we adopt ANSI C as our standard. You might wantto look ahead to §1.2 where ANSI C function prototypes are discussed in more detail.Note our convention of handling all errors and exceptional cases with a statementlike nrerror("some error message");. The function nrerror() is part of asmall file of utility programs, nrutil.c, listed in Appendix B at the back of thebook.

This Appendix includes a number of other utilities that we will describe later inthis chapter. Function nrerror() prints the indicated error message to your stderrdevice (usually your terminal screen), and then invokes the function exit(), whichterminates execution. The function exit() is in every C library we know of; but ifyou find it missing, you can modify nrerror() so that it does anything else that willhalt execution. For example, you can have it pause for input from the keyboard, andthen manually interrupt execution.

In some applications, you will want to modifynrerror() to do more sophisticated error handling, for example to transfer controlsomewhere else, with an error flag or error code set.We will have more to say about the C programming language, its conventionsand style, in §1.1 and §1.2.Computational Environment and Program ValidationOur goal is that the programs in this book be as portable as possible, acrossdifferent platforms (models of computer), across different operating systems, andacross different C compilers. C was designed with this type of portability inmind.

Nevertheless, we have found that there is no substitute for actually checkingall programs on a variety of compilers, in the process uncovering differences inlibrary structure or contents, and even occasional differences in allowed syntax. Assurrogates for the large number of possible combinations, we have tested all theprograms in this book on the combinations of machines, operating systems, andcompilers shown on the accompanying table. More generally, the programs shouldrun without modification on any compiler that implements the ANSI C standard,as described for example in Harbison and Steele’s excellent book [1].

With smallmodifications, our programs should run on any compiler that implements the older,de facto K&R standard [2]. An example of the kind of trivial incompatibility towatch out for is that ANSI C requires the memory allocation functions malloc()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).t=c/1236.85;t2=t*t;as=359.2242+29.105356*c;You aren’t really intended to understandam=306.0253+385.816918*c+0.010730*t2;this algorithm, but it does work!*jd=2415020+28L*n+7L*nph;xtra=0.75933+1.53058868*c+((1.178e-4)-(1.55e-7)*t)*t2;if (nph == 0 || nph == 2)xtra += (0.1734-3.93e-4*t)*sin(RAD*as)-0.4068*sin(RAD*am);else if (nph == 1 || nph == 3)xtra += (0.1721-4.0e-4*t)*sin(RAD*as)-0.6280*sin(RAD*am);else nrerror("nph is unknown in flmoon");This is how we will indicate errori=(int)(xtra >= 0.0 ? floor(xtra) : ceil(xtra-1.0));conditions.*jd += i;*frac=xtra-i;31.0 IntroductionTested Machines and CompilersHardwareO/S VersionCompiler VersionMS-DOS 5.0/Windows 3.1Microsoft C/C++ 7.0IBM PC compatible 486/33MS-DOS 5.0Borland C/C++ 2.0IBM RS/6000AIX 3.2IBM xlc 1.02DECstation 5000/25ULTRIX 4.2ACodeCenter (Saber) C 3.1.1DECsystem 5400ULTRIX 4.1GNU C Compiler 2.1Sun SPARCstation 2SunOS 4.1GNU C Compiler 1.40DECstation 5000/200ULTRIX 4.2DEC RISC C 2.1*Sun SPARCstation 2SunOS 4.1Sun cc 1.1**compiler version does not fully implement ANSI C; only K&R validatedand free() to be declared via the header stdlib.h; some older compilers requirethem to be declared with the header file malloc.h, while others regard them asinherent in the language and require no header file at all.In validating the programs, we have taken the program source code directlyfrom the machine-readable form of the book’s manuscript, to decrease the chanceof propagating typographical errors.

“Driver” or demonstration programs that weused as part of our validations are available separately as the Numerical RecipesExample Book (C), as well as in machine-readable form. If you plan to use morethan a few of the programs in this book, or if you plan to use programs in this bookon more than one different computer, then you may find it useful to obtain a copyof these demonstration programs.Of course we would be foolish to claim that there are no bugs in our programs,and we do not make such a claim. We have been very careful, and have benefittedfrom the experience of the many readers who have written to us.

If you find a newbug, please document it and tell us!Compatibility with the First EditionIf you are accustomed to the Numerical Recipes routines of the First Edition, restassured: almost all of them are still here, with the same names and functionalities,often with major improvements in the code itself. In addition, we hope that youwill soon become equally familiar with the added capabilities of the more than 100routines that are new to this edition.We have retired a small number of First Edition routines, those that we believeto be clearly dominated by better methods implemented in this edition.

A table,following, lists the retired routines and suggests replacements.First Edition users should also be aware that some routines common toboth editions have alterations in their calling interfaces, so are not directly “plugcompatible.” A fairly complete list is: chsone, chstwo, covsrt, dfpmin, laguer,lfit, memcof, mrqcof, mrqmin, pzextr, ran4, realft, rzextr, shoot, shootf.There may be others (depending in part on which printing of the First Edition is takenfor the comparison).