Numerical Recipes in C

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 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).

Numerical Recipes in C
The Art of Scientiﬁc Computing

Cambridge New York Port Chester Melbourne Sydney EXXON Research and Engineering Company Harvard-Smithsonian Center for Astrophysics Department of Physics, Cornell University

CAMBRIDGE UNIVERSITY PRESS

William T. Vetterling Saul A. Teukolsky

Brian P. Flannery

Second Edition

William H. Press

Polaroid Corporation

Published by the Press Syndicate of the University of Cambridge The Pitt Building, Trumpington Street, Cambridge CB2 1RP 40 West 20th Street, New York, NY 10011-4211, USA 477 Williamstown Road, Port Melbourne, VIC, 3207, Australia Copyright c Cambridge University Press 1988, 1992 except for §13.10 and Appendix B, which are placed into the public domain, and except for all other computer programs and procedures, which are Copyright c Numerical Recipes Software 1987, 1988, 1992, 1997, 2002 All Rights Reserved. Some sections of this book were originally published, in different form, in Computers in Physics magazine, Copyright c American Institute of Physics, 1988–1992. First Edition originally published 1988; Second Edition originally published 1992. Reprinted with corrections, 1993, 1994, 1995, 1997, 2002. This reprinting is corrected to software version 2.10 Printed in the United States of America Typeset in TEX

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 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).

Without an additional license to use the contained software, this book is intended as a text and reference book, for reading purposes only. A free license for limited use of the software by the individual owner of a copy of this book who personally types one or more routines into a single computer is granted under terms described on p. xvii. See the section “License Information” (pp. xvi–xviii) for information on obtaining more general licenses at low cost. Machine-readable media containing the software in this book, with included licenses for use on a single screen, are available from Cambridge University Press. See the order form at the back of the book, email to “orders@cup.org” (North America) or “directcustserv@cambridge.org” (rest of world), or write to Cambridge University Press, 110 Midland Avenue, Port Chester, NY 10573 (USA), for further information. The software may also be downloaded, with immediate purchase of a license also possible, from the Numerical Recipes Software Web Site (http://www.nr.com). Unlicensed transfer of Numerical Recipes programs to any other format, or to any computer except one that is speciﬁcally licensed, is strictly prohibited. Technical questions, corrections, and requests for information should be addressed to Numerical Recipes Software, P.O. Box 380243, Cambridge, MA 02238-0243 (USA), email “info@nr.com”, or fax 781 863-1739.

Library of Congress Cataloging in Publication Data Numerical recipes in C : the art of scientiﬁc computing / William H. Press . . . [et al.]. – 2nd ed. Includes bibliographical references (p. ) and index. ISBN 0-521-43108-5 1. Numerical analysis–Computer programs. 2. Science–Mathematics–Computer programs. 3. C (Computer program language) I. Press, William H. QA297.N866 1992 92-8876 519.4 0285 53–dc20 A catalog record for this book is available from the British Library. ISBN ISBN ISBN ISBN ISBN 0 0 0 0 0 521 521 521 521 521 43108 43720 75037 75035 75036 5 2 7 0 9 Book Example book in C C/C++ CDROM (Windows/Macintosh) Complete CDROM (Windows/Macintosh) Complete CDROM (UNIX/Linux)

Contents

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 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).

Preface to the Second Edition Preface to the First Edition License Information Computer Programs by Chapter and Section 1 Preliminaries
1.0 Introduction 1.1 Program Organization and Control Structures 1.2 Some C Conventions for Scientiﬁc Computing 1.3 Error, Accuracy, and Stability

xi xiv xvi xix 1
1 5 15 28

2

Solution of Linear Algebraic Equations
2.0 Introduction 2.1 Gauss-Jordan Elimination 2.2 Gaussian Elimination with Backsubstitution 2.3 LU Decomposition and Its Applications 2.4 Tridiagonal and Band Diagonal Systems of Equations 2.5 Iterative Improvement of a Solution to Linear Equations 2.6 Singular Value Decomposition 2.7 Sparse Linear Systems 2.8 Vandermonde Matrices and Toeplitz Matrices 2.9 Cholesky Decomposition 2.10 QR Decomposition 2.11 Is Matrix Inversion an N 3 Process?

32
32 36 41 43 50 55 59 71 90 96 98 102

3

Interpolation and Extrapolation
3.0 Introduction 3.1 Polynomial Interpolation and Extrapolation 3.2 Rational Function Interpolation and Extrapolation 3.3 Cubic Spline Interpolation 3.4 How to Search an Ordered Table 3.5 Coefﬁcients of the Interpolating Polynomial 3.6 Interpolation in Two or More Dimensions v

105
105 108 111 113 117 120 123

vi

Contents

4

Integration of Functions
4.0 Introduction 4.1 Classical Formulas for Equally Spaced Abscissas 4.2 Elementary Algorithms 4.3 Romberg Integration 4.4 Improper Integrals 4.5 Gaussian Quadratures and Orthogonal Polynomials 4.6 Multidimensional Integrals

129
129 130 136 140 141 147 161

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 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).

5

Evaluation of Functions
5.0 Introduction 5.1 Series and Their Convergence 5.2 Evaluation of Continued Fractions 5.3 Polynomials and Rational Functions 5.4 Complex Arithmetic 5.5 Recurrence Relations and Clenshaw’s Recurrence Formula 5.6 Quadratic and Cubic Equations 5.7 Numerical Derivatives 5.8 Chebyshev Approximation 5.9 Derivatives or Integrals of a Chebyshev-approximated Function 5.10 Polynomial Approximation from Chebyshev Coefﬁcients 5.11 Economization of Power Series 5.12 Pad´ Approximants e 5.13 Rational Chebyshev Approximation 5.14 Evaluation of Functions by Path Integration

165
165 165 169 173 176 178 183 186 190 195 197 198 200 204 208

6

Special Functions
6.0 Introduction 6.1 Gamma Function, Beta Function, Factorials, Binomial Coefﬁcients 6.2 Incomplete Gamma Function, Error Function, Chi-Square Probability Function, Cumulative Poisson Function 6.3 Exponential Integrals 6.4 Incomplete Beta Function, Student’s Distribution, F-Distribution, Cumulative Binomial Distribution 6.5 Bessel Functions of Integer Order 6.6 Modiﬁed Bessel Functions of Integer Order 6.7 Bessel Functions of Fractional Order, Airy Functions, Spherical Bessel Functions 6.8 Spherical Harmonics 6.9 Fresnel Integrals, Cosine and Sine Integrals 6.10 Dawson’s Integral 6.11 Elliptic Integrals and Jacobian Elliptic Functions 6.12 Hypergeometric Functions

212
212 213 216 222 226 230 236 240 252 255 259 261 271

7

Random Numbers
7.0 Introduction 7.1 Uniform Deviates

274
274 275

Contents

vii 287 290 296 300 304 309 316

7.2 Transformation Method: Exponential and Normal Deviates 7.3 Rejection Method: Gamma, Poisson, Binomial Deviates 7.4 Generation of Random Bits 7.5 Random Sequences Based on Data Encryption 7.6 Simple Monte Carlo Integration 7.7 Quasi- (that is, Sub-) Random Sequences 7.8 Adaptive and Recursive Monte Carlo Methods

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 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).

8

Sorting
8.0 Introduction 8.1 Straight Insertion and Shell’s Method 8.2 Quicksort 8.3 Heapsort 8.4 Indexing and Ranking 8.5 Selecting the M th Largest 8.6 Determination of Equivalence Classes

329
329 330 332 336 338 341 345

9

Root Finding and Nonlinear Sets of Equations
9.0 Introduction 9.1 Bracketing and Bisection 9.2 Secant Method, False Position Method, and Ridders’ Method 9.3 Van Wijngaarden–Dekker–Brent Method 9.4 Newton-Raphson Method Using Derivative 9.5 Roots of Polynomials 9.6 Newton-Raphson Method for Nonlinear Systems of Equations 9.7 Globally Convergent Methods for Nonlinear Systems of Equations

347
347 350 354 359 362 369 379 383

10 Minimization or Maximization of Functions
10.0 Introduction 10.1 Golden Section Search in One Dimension 10.2 Parabolic Interpolation and Brent’s Method in One Dimension 10.3 One-Dimensional Search with First Derivatives 10.4 Downhill Simplex Method in Multidimensions 10.5 Direction Set (Powell’s) Methods in Multidimensions 10.6 Conjugate Gradient Methods in Multidimensions 10.7 Variable Metric Methods in Multidimensions 10.8 Linear Programming and the Simplex Method 10.9 Simulated Annealing Methods

394
394 397 402 405 408 412 420 425 430 444

11 Eigensystems
11.0 Introduction 11.1 Jacobi Transformations of a Symmetric Matrix 11.2 Reduction of a Symmetric Matrix to Tridiagonal Form: Givens and Householder Reductions 11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 11.4 Hermitian Matrices 11.5 Reduction of a General Matrix to Hessenberg Form

456
456 463 469 475 481 482

viii

Contents

11.6 The QR Algorithm for Real Hessenberg Matrices 11.7 Improving Eigenvalues and/or Finding Eigenvectors by Inverse Iteration

486 493

12 Fast Fourier Transform
12.0 Introduction 12.1 Fourier Transform of Discretely Sampled Data 12.2 Fast Fourier Transform (FFT) 12.3 FFT of Real Functions, Sine and Cosine Transforms 12.4 FFT in Two or More Dimensions 12.5 Fourier Transforms of Real Data in Two and Three Dimensions 12.6 External Storage or Memory-Local FFTs

496
496 500 504 510 521 525 532
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 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).

13 Fourier and Spectral Applications
13.0 Introduction 13.1 Convolution and Deconvolution Using the FFT 13.2 Correlation and Autocorrelation Using the FFT 13.3 Optimal (Wiener) Filtering with the FFT 13.4 Power Spectrum Estimation Using the FFT 13.5 Digital Filtering in the Time Domain 13.6 Linear Prediction and Linear Predictive Coding 13.7 Power Spectrum Estimation by the Maximum Entropy (All Poles) Method 13.8 Spectral Analysis of Unevenly Sampled Data 13.9 Computing Fourier Integrals Using the FFT 13.10 Wavelet Transforms 13.11 Numerical Use of the Sampling Theorem

537
537 538 545 547 549 558 564 572 575 584 591 606

14 Statistical Description of Data
14.0 Introduction 14.1 Moments of a Distribution: Mean, Variance, Skewness, and So Forth 14.2 Do Two Distributions Have the Same Means or Variances? 14.3 Are Two Distributions Different? 14.4 Contingency Table Analysis of Two Distributions 14.5 Linear Correlation 14.6 Nonparametric or Rank Correlation 14.7 Do Two-Dimensional Distributions Differ? 14.8 Savitzky-Golay Smoothing Filters

609
609 610 615 620 628 636 639 645 650

15 Modeling of Data
15.0 Introduction 15.1 Least Squares as a Maximum Likelihood Estimator 15.2 Fitting Data to a Straight Line 15.3 Straight-Line Data with Errors in Both Coordinates 15.4 General Linear Least Squares 15.5 Nonlinear Models

656
656 657 661 666 671 681

Contents

ix 689 699

15.6 Conﬁdence Limits on Estimated Model Parameters 15.7 Robust Estimation

16 Integration of Ordinary Differential Equations
16.0 Introduction 16.1 Runge-Kutta Method 16.2 Adaptive Stepsize Control for Runge-Kutta 16.3 Modiﬁed Midpoint Method 16.4 Richardson Extrapolation and the Bulirsch-Stoer Method 16.5 Second-Order Conservative Equations 16.6 Stiff Sets of Equations 16.7 Multistep, Multivalue, and Predictor-Corrector Methods

707
707 710 714 722 724 732 734 747
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 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).

17 Two Point Boundary Value Problems
17.0 Introduction 17.1 The Shooting Method 17.2 Shooting to a Fitting Point 17.3 Relaxation Methods 17.4 A Worked Example: Spheroidal Harmonics 17.5 Automated Allocation of Mesh Points 17.6 Handling Internal Boundary Conditions or Singular Points

753
753 757 760 762 772 783 784

18 Integral Equations and Inverse Theory
18.0 Introduction 18.1 Fredholm Equations of the Second Kind 18.2 Volterra Equations 18.3 Integral Equations with Singular Kernels 18.4 Inverse Problems and the Use of A Priori Information 18.5 Linear Regularization Methods 18.6 Backus-Gilbert Method 18.7 Maximum Entropy Image Restoration

788
788 791 794 797 804 808 815 818

19 Partial Differential Equations
19.0 Introduction 19.1 Flux-Conservative Initial Value Problems 19.2 Diffusive Initial Value Problems 19.3 Initial Value Problems in Multidimensions 19.4 Fourier and Cyclic Reduction Methods for Boundary Value Problems 19.5 Relaxation Methods for Boundary Value Problems 19.6 Multigrid Methods for Boundary Value Problems

827
827 834 847 853 857 863 871

20 Less-Numerical Algorithms
20.0 Introduction 20.1 Diagnosing Machine Parameters 20.2 Gray Codes

889
889 889 894

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 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).

926

930

951 Index of Programs and Dependencies General Index 965

940 Appendix C: Complex Arithmetic Appendix B: Utility Routines 948

896 903 910 915

x

Appendix A: Table of Prototype Declarations
20.3 Cyclic Redundancy and Other Checksums 20.4 Huffman Coding and Compression of Data 20.5 Arithmetic Coding 20.6 Arithmetic at Arbitrary Precision

Contents

References

Preface to the Second Edition
Our aim in writing the original edition of Numerical Recipes was to provide a book that combined general discussion, analytical mathematics, algorithmics, and actual working programs. The success of the ﬁrst edition puts us now in a difﬁcult, though hardly unenviable, position. We wanted, then and now, to write a book that is informal, fearlessly editorial, unesoteric, and above all useful. There is a danger that, if we are not careful, we might produce a second edition that is weighty, balanced, scholarly, and boring. It is a mixed blessing that we know more now than we did six years ago. Then, we were making educated guesses, based on existing literature and our own research, about which numerical techniques were the most important and robust. Now, we have the beneﬁt of direct feedback from a large reader community. Letters to our alter-ego enterprise, Numerical Recipes Software, are in the thousands per year. (Please, don’t telephone us.) Our post ofﬁce box has become a magnet for letters pointing out that we have omitted some particular technique, well known to be important in a particular ﬁeld of science or engineering. We value such letters, and digest them carefully, especially when they point us to speciﬁc references in the literature. The inevitable result of this input is that this Second Edition of Numerical Recipes is substantially larger than its predecessor, in fact about 50% larger both in words and number of included programs (the latter now numbering well over 300). “Don’t let the book grow in size,” is the advice that we received from several wise colleagues. We have tried to follow the intended spirit of that advice, even as we violate the letter of it. We have not lengthened, or increased in difﬁculty, the book’s principal discussions of mainstream topics. Many new topics are presented at this same accessible level. Some topics, both from the earlier edition and new to this one, are now set in smaller type that labels them as being “advanced.” The reader who ignores such advanced sections completely will not, we think, ﬁnd any lack of continuity in the shorter volume that results. Here are some highlights of the new material in this Second Edition: • a new chapter on integral equations and inverse methods • a detailed treatment of multigrid methods for solving elliptic partial differential equations • routines for band diagonal linear systems • improved routines for linear algebra on sparse matrices • Cholesky and QR decomposition • orthogonal polynomials and Gaussian quadratures for arbitrary weight functions • methods for calculating numerical derivatives • Pad´ approximants, and rational Chebyshev approximation e • Bessel functions, and modiﬁed Bessel functions, of fractional order; and several other new special functions • improved random number routines • quasi-random sequences • routines for adaptive and recursive Monte Carlo integration in highdimensional spaces • globally convergent methods for sets of nonlinear equations xi

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 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).

xii • • • • • • • • • • • • • •

Preface to the Second Edition

simulated annealing minimization for continuous control spaces fast Fourier transform (FFT) for real data in two and three dimensions fast Fourier transform (FFT) using external storage improved fast cosine transform routines wavelet transforms Fourier integrals with upper and lower limits spectral analysis on unevenly sampled data Savitzky-Golay smoothing ﬁlters ﬁtting straight line data with errors in both coordinates a two-dimensional Kolmogorov-Smirnoff test the statistical bootstrap method embedded Runge-Kutta-Fehlberg methods for differential equations high-order methods for stiff differential equations a new chapter on “less-numerical” algorithms, including Huffman and arithmetic coding, arbitrary precision arithmetic, and several other topics. Consult the Preface to the First Edition, following, or the Table of Contents, for a list of the more “basic” subjects treated.

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 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).

Acknowledgments
It is not possible for us to list by name here all the readers who have made useful suggestions; we are grateful for these. In the text, we attempt to give speciﬁc attribution for ideas that appear to be original, and not known in the literature. We apologize in advance for any omissions. Some readers and colleagues have been particularly generous in providing us with ideas, comments, suggestions, and programs for this Second Edition. We especially want to thank George Rybicki, Philip Pinto, Peter Lepage, Robert Lupton, Douglas Eardley, Ramesh Narayan, David Spergel, Alan Oppenheim, Sallie Baliunas, Scott Tremaine, Glennys Farrar, Steven Block, John Peacock, Thomas Loredo, Matthew Choptuik, Gregory Cook, L. Samuel Finn, P. Deuﬂhard, Harold Lewis, Peter Weinberger, David Syer, Richard Ferch, Steven Ebstein, Bradley Keister, and William Gould. We have been helped by Nancy Lee Snyder’s mastery of a complicated TEX manuscript. We express appreciation to our editors Lauren Cowles and Alan Harvey at Cambridge University Press, and to our production editor Russell Hahn. We remain, of course, grateful to the individuals acknowledged in the Preface to the First Edition. Special acknowledgment is due to programming consultant Seth Finkelstein, who wrote, rewrote, or inﬂuenced many of the routines in this book, as well as in its FORTRAN-language twin and the companion Example books. Our project has beneﬁted enormously from Seth’s talent for detecting, and following the trail of, even very slight anomalies (often compiler bugs, but occasionally our errors), and from his good programming sense. To the extent that this edition of Numerical Recipes in C has a more graceful and “C-like” programming style than its predecessor, most of the credit goes to Seth. (Of course, we accept the blame for the FORTRANish lapses that still remain.) We prepared this book for publication on DEC and Sun workstations running the UNIX operating system, and on a 486/33 PC compatible running MS-DOS 5.0/Windows 3.0. (See §1.0 for a list of additional computers used in

Preface to the Second Edition

xiii

program tests.) We enthusiastically recommend the principal software used: GNU Emacs, TEX, Perl, Adobe Illustrator, and PostScript. Also used were a variety of C compilers – too numerous (and sometimes too buggy) for individual acknowledgment. It is a sobering fact that our standard test suite (exercising all the routines in this book) has uncovered compiler bugs in many of the compilers tried. When possible, we work with developers to see that such bugs get ﬁxed; we encourage interested compiler developers to contact us about such arrangements. WHP and SAT acknowledge the continued support of the U.S. National Science Foundation for their research on computational methods. D.A.R.P.A. support is acknowledged for §13.10 on wavelets. June, 1992 William H. Press Saul A. Teukolsky William T. Vetterling Brian P. Flannery

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 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).

Preface to the First Edition
We call this book Numerical Recipes for several reasons. In one sense, this book is indeed a “cookbook” on numerical computation. However there is an important distinction between a cookbook and a restaurant menu. The latter presents choices among complete dishes in each of which the individual ﬂavors are blended and disguised. The former — and this book — reveals the individual ingredients and explains how they are prepared and combined. Another purpose of the title is to connote an eclectic mixture of presentational techniques. This book is unique, we think, in offering, for each topic considered, a certain amount of general discussion, a certain amount of analytical mathematics, a certain amount of discussion of algorithmics, and (most important) actual implementations of these ideas in the form of working computer routines. Our task has been to ﬁnd the right balance among these ingredients for each topic. You will ﬁnd that for some topics we have tilted quite far to the analytic side; this where we have felt there to be gaps in the “standard” mathematical training. For other topics, where the mathematical prerequisites are universally held, we have tilted towards more in-depth discussion of the nature of the computational algorithms, or towards practical questions of implementation. We admit, therefore, to some unevenness in the “level” of this book. About half of it is suitable for an advanced undergraduate course on numerical computation for science or engineering majors. The other half ranges from the level of a graduate course to that of a professional reference. Most cookbooks have, after all, recipes at varying levels of complexity. An attractive feature of this approach, we think, is that the reader can use the book at increasing levels of sophistication as his/her experience grows. Even inexperienced readers should be able to use our most advanced routines as black boxes. Having done so, we hope that these readers will subsequently go back and learn what secrets are inside. If there is a single dominant theme in this book, it is that practical methods of numerical computation can be simultaneously efﬁcient, clever, and — important — clear. The alternative viewpoint, that efﬁcient computational methods must necessarily be so arcane and complex as to be useful only in “black box” form, we ﬁrmly reject. Our purpose in this book is thus to open up a large number of computational black boxes to your scrutiny. We want to teach you to take apart these black boxes and to put them back together again, modifying them to suit your speciﬁc needs. We assume that you are mathematically literate, i.e., that you have the normal mathematical preparation associated with an undergraduate degree in a physical science, or engineering, or economics, or a quantitative social science. We assume that you know how to program a computer. We do not assume that you have any prior formal knowledge of numerical analysis or numerical methods. The scope of Numerical Recipes is supposed to be “everything up to, but not including, partial differential equations.” We honor this in the breach: First, we do have one introductory chapter on methods for partial differential equations (Chapter 19). Second, we obviously cannot include everything else. All the so-called “standard” topics of a numerical analysis course have been included in this book: xiv

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 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).

Preface to the First Edition

xv

linear equations (Chapter 2), interpolation and extrapolation (Chaper 3), integration (Chaper 4), nonlinear root-ﬁnding (Chapter 9), eigensystems (Chapter 11), and ordinary differential equations (Chapter 16). Most of these topics have been taken beyond their standard treatments into some advanced material which we have felt to be particularly important or useful. Some other subjects that we cover in detail are not usually found in the standard numerical analysis texts. These include the evaluation of functions and of particular special functions of higher mathematics (Chapters 5 and 6); random numbers and Monte Carlo methods (Chapter 7); sorting (Chapter 8); optimization, including multidimensional methods (Chapter 10); Fourier transform methods, including FFT methods and other spectral methods (Chapters 12 and 13); two chapters on the statistical description and modeling of data (Chapters 14 and 15); and two-point boundary value problems, both shooting and relaxation methods (Chapter 17). The programs in this book are included in ANSI-standard C. Versions of the book in FORTRAN, Pascal, and BASIC are available separately. We have more to say about the C language, and the computational environment assumed by our routines, in §1.1 (Introduction).

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 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).

Acknowledgments
Many colleagues have been generous in giving us the beneﬁt of their numerical and computational experience, in providing us with programs, in commenting on the manuscript, or in general encouragement. We particularly wish to thank George Rybicki, Douglas Eardley, Philip Marcus, Stuart Shapiro, Paul Horowitz, Bruce Musicus, Irwin Shapiro, Stephen Wolfram, Henry Abarbanel, Larry Smarr, Richard Muller, John Bahcall, and A.G.W. Cameron. We also wish to acknowledge two individuals whom we have never met: Forman Acton, whose 1970 textbook Numerical Methods that Work (New York: Harper and Row) has surely left its stylistic mark on us; and Donald Knuth, both for his series of books on The Art of Computer Programming (Reading, MA: Addison-Wesley), and for TEX, the computer typesetting language which immensely aided production of this book. Research by the authors on computational methods was supported in part by the U.S. National Science Foundation. October, 1985 William H. Press Brian P. Flannery Saul A. Teukolsky William T. Vetterling

License Information
Read this section if you want to use the programs in this book on a computer. You’ll need to read the following Disclaimer of Warranty, get the programs onto your computer, and acquire a Numerical Recipes software license. (Without this license, which can be the free “immediate license” under terms described below, the book is intended as a text and reference book, for reading purposes only.)

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 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).

Disclaimer of Warranty
We make no warranties, express or implied, that the programs contained in this volume are free of error, or are consistent with any particular standard of merchantability, or that they will meet your requirements for any particular application. They should not be relied on for solving a problem whose incorrect solution could result in injury to a person or loss of property. If you do use the programs in such a manner, it is at your own risk. The authors and publisher disclaim all liability for direct or consequential damages resulting from your use of the programs.

How to Get the Code onto Your Computer
Pick one of the following methods:
• You can type the programs from this book directly into your computer. In this case, the only kind of license available to you is the free “immediate license” (see below). You are not authorized to transfer or distribute a machine-readable copy to any other person, nor to have any other person type the programs into a computer on your behalf. We do not want to hear bug reports from you if you choose this option, because experience has shown that virtually all reported bugs in such cases are typing errors! • You can download the Numerical Recipes programs electronically from the Numerical Recipes On-Line Software Store, located at http://www.nr.com, our Web site. All the ﬁles (Recipes and demonstration programs) are packaged as a single compressed ﬁle. You’ll need to purchase a license to download and unpack them. Any number of single-screen licenses can be purchased instantly (with discount for multiple screens) from the On-Line Store, with fees that depend on your operating system (Windows or Macintosh versus Linux or UNIX) and whether you are afﬁliated with an educational institution. Purchasing a singlescreen license is also the way to start if you want to acquire a more general (site or corporate) license; your single-screen cost will be subtracted from the cost of any later license upgrade. • You can purchase media containing the programs from Cambridge University Press. A CD-ROM version in ISO-9660 format for Windows and Macintosh systems contains the complete C software, and also the C++ version. More extensive CDROMs in ISO-9660 format for Windows, Macintosh, and UNIX/Linux systems are also available; these include the C, C++, and Fortran versions on a single CD-ROM (as well as versions in Pascal and BASIC from the ﬁrst edition). These CD-ROMs are available with a single-screen license for Windows or Macintosh (order ISBN 0 521 750350), or (at a slightly higher price) with a single-screen license for UNIX/Linux workstations (order ISBN 0 521 750369). Orders for media from Cambridge University Press can be placed at 800 872-7423 (North America only) or by email to orders@cup.org (North America) or directcustserv@cambridge.org (rest of world). Or, visit the Web site http://www.cambridge.org.

xvi

License Information

xvii

Types of License Offered
Here are the types of licenses that we offer. Note that some types are automatically acquired with the purchase of media from Cambridge University Press, or of an unlocking password from the Numerical Recipes On-Line Software Store, while other types of licenses require that you communicate speciﬁcally with Numerical Recipes Software (email: orders@nr.com or fax: 781 863-1739). Our Web site http://www.nr.com has additional information.
• [“Immediate License”] If you are the individual owner of a copy of this book and you type one or more of its routines into your computer, we authorize you to use them on that computer for your own personal and noncommercial purposes. You are not authorized to transfer or distribute machine-readable copies to any other person, or to use the routines on more than one machine, or to distribute executable programs containing our routines. This is the only free license. • [“Single-Screen License”] This is the most common type of low-cost license, with terms governed by our Single Screen (Shrinkwrap) License document (complete terms available through our Web site). Basically, this license lets you use Numerical Recipes routines on any one screen (PC, workstation, X-terminal, etc.). You may also, under this license, transfer pre-compiled, executable programs incorporating our routines to other, unlicensed, screens or computers, providing that (i) your application is noncommercial (i.e., does not involve the selling of your program for a fee), (ii) the programs were ﬁrst developed, compiled, and successfully run on a licensed screen, and (iii) our routines are bound into the programs in such a manner that they cannot be accessed as individual routines and cannot practicably be unbound and used in other programs. That is, under this license, your program user must not be able to use our programs as part of a program library or “mix-andmatch” workbench. Conditions for other types of commercial or noncommercial distribution may be found on our Web site (http://www.nr.com). • [“Multi-Screen, Server, Site, and Corporate Licenses”] The terms of the Single Screen License can be extended to designated groups of machines, deﬁned by number of screens, number of machines, locations, or ownership. Signiﬁcant discounts from the corresponding single-screen prices are available when the estimated number of screens exceeds 40. Contact Numerical Recipes Software (email: orders@nr.com or fax: 781 863-1739) for details. • [“Course Right-to-Copy License”] Instructors at accredited educational institutions who have adopted this book for a course, and who have already purchased a Single Screen License (either acquired with the purchase of media, or from the Numerical Recipes On-Line Software Store), may license the programs for use in that course as follows: Mail your name, title, and address; the course name, number, dates, and estimated enrollment; and advance payment of $5 per (estimated) student to Numerical Recipes Software, at this address: P.O. Box 243, Cambridge, MA 02238 (USA). You will receive by return mail a license authorizing you to make copies of the programs for use by your students, and/or to transfer the programs to a machine accessible to your students (but only for the duration of the course).

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 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).

About Copyrights on Computer Programs
Like artistic or literary compositions, computer programs are protected by copyright. Generally it is an infringement for you to copy into your computer a program from a copyrighted source. (It is also not a friendly thing to do, since it deprives the program’s author of compensation for his or her creative effort.) Under

xviii

License Information

copyright law, all “derivative works” (modiﬁed versions, or translations into another computer language) also come under the same copyright as the original work. Copyright does not protect ideas, but only the expression of those ideas in a particular form. In the case of a computer program, the ideas consist of the program’s methodology and algorithm, including the necessary sequence of steps adopted by the programmer. The expression of those ideas is the program source code (particularly any arbitrary or stylistic choices embodied in it), its derived object code, and any other derivative works. If you analyze the ideas contained in a program, and then express those ideas in your own completely different implementation, then that new program implementation belongs to you. That is what we have done for those programs in this book that are not entirely of our own devising. When programs in this book are said to be “based” on programs published in copyright sources, we mean that the ideas are the same. The expression of these ideas as source code is our own. We believe that no material in this book infringes on an existing copyright.

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 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).

Trademarks
Several registered trademarks appear within the text of this book: Sun is a trademark of Sun Microsystems, Inc. SPARC and SPARCstation are trademarks of SPARC International, Inc. Microsoft, Windows 95, Windows NT, PowerStation, and MS are trademarks of Microsoft Corporation. DEC, VMS, Alpha AXP, and ULTRIX are trademarks of Digital Equipment Corporation. IBM is a trademark of International Business Machines Corporation. Apple and Macintosh are trademarks of Apple Computer, Inc. UNIX is a trademark licensed exclusively through X/Open Co. Ltd. IMSL is a trademark of Visual Numerics, Inc. NAG refers to proprietary computer software of Numerical Algorithms Group (USA) Inc. PostScript and Adobe Illustrator are trademarks of Adobe Systems Incorporated. Last, and no doubt least, Numerical Recipes (when identifying products) is a trademark of Numerical Recipes Software.

Attributions
The fact that ideas are legally “free as air” in no way supersedes the ethical requirement that ideas be credited to their known originators. When programs in this book are based on known sources, whether copyrighted or in the public domain, published or “handed-down,” we have attempted to give proper attribution. Unfortunately, the lineage of many programs in common circulation is often unclear. We would be grateful to readers for new or corrected information regarding attributions, which we will attempt to incorporate in subsequent printings.

Computer Programs by Chapter and Section
1.0 1.1 1.1 1.1 2.1 2.3 2.3 2.4 2.4 2.4 2.4 2.5 2.6 2.6 2.6 2.7 2.7 2.7 2.7 2.7 2.7 2.7 2.7 2.7 2.7 2.7 2.8 2.8 2.9 2.9 2.10 2.10 2.10 2.10 2.10 3.1 3.2 3.3 3.3 3.4 flmoon julday badluk caldat gaussj ludcmp lubksb tridag banmul bandec banbks mprove svbksb svdcmp pythag cyclic sprsin sprsax sprstx sprstp sprspm sprstm linbcg snrm atimes asolve vander toeplz choldc cholsl qrdcmp qrsolv rsolv qrupdt rotate polint ratint spline splint locate calculate phases of the moon by date Julian Day number from calendar date Friday the 13th when the moon is full calendar date from Julian day number Gauss-Jordan matrix inversion and linear equation solution linear equation solution, LU decomposition linear equation solution, backsubstitution solution of tridiagonal systems multiply vector by band diagonal matrix band diagonal systems, decomposition band diagonal systems, backsubstitution linear equation solution, iterative improvement singular value backsubstitution singular value decomposition of a matrix calculate (a2 + b2 )1/2 without overﬂow solution of cyclic tridiagonal systems convert matrix to sparse format product of sparse matrix and vector product of transpose sparse matrix and vector transpose of sparse matrix pattern multiply two sparse matrices threshold multiply two sparse matrices biconjugate gradient solution of sparse systems used by linbcg for vector norm used by linbcg for sparse multiplication used by linbcg for preconditioner solve Vandermonde systems solve Toeplitz systems Cholesky decomposition Cholesky backsubstitution QR decomposition QR backsubstitution right triangular backsubstitution update a QR decomposition Jacobi rotation used by qrupdt polynomial interpolation rational function interpolation construct a cubic spline cubic spline interpolation search an ordered table by bisection xix
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 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).

xx 3.4 3.5 3.5 3.6 3.6 3.6 3.6 3.6 4.2 4.2 4.2 4.3 4.4 4.4 4.4 4.4 4.4 4.4 4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.6 5.1 5.3 5.3 5.3 5.7 5.8 5.8 5.9 5.9 5.10 5.10 5.11 5.12 5.13 6.1 6.1 6.1 6.1 hunt polcoe polcof polin2 bcucof bcuint splie2 splin2 trapzd qtrap qsimp qromb midpnt qromo midinf midsql midsqu midexp qgaus gauleg gaulag gauher gaujac gaucof orthog quad3d eulsum ddpoly poldiv ratval dfridr chebft chebev chder chint chebpc pcshft pccheb pade ratlsq gammln factrl bico factln

Computer Programs by Chapter and Section

search a table when calls are correlated polynomial coefﬁcients from table of values polynomial coefﬁcients from table of values two-dimensional polynomial interpolation construct two-dimensional bicubic two-dimensional bicubic interpolation construct two-dimensional spline two-dimensional spline interpolation trapezoidal rule integrate using trapezoidal rule integrate using Simpson’s rule integrate using Romberg adaptive method extended midpoint rule integrate using open Romberg adaptive method integrate a function on a semi-inﬁnite interval integrate a function with lower square-root singularity integrate a function with upper square-root singularity integrate a function that decreases exponentially integrate a function by Gaussian quadratures Gauss-Legendre weights and abscissas Gauss-Laguerre weights and abscissas Gauss-Hermite weights and abscissas Gauss-Jacobi weights and abscissas quadrature weights from orthogonal polynomials construct nonclassical orthogonal polynomials integrate a function over a three-dimensional space sum a series by Euler–van Wijngaarden algorithm evaluate a polynomial and its derivatives divide one polynomial by another evaluate a rational function numerical derivative by Ridders’ method ﬁt a Chebyshev polynomial to a function Chebyshev polynomial evaluation derivative of a function already Chebyshev ﬁtted integrate a function already Chebyshev ﬁtted polynomial coefﬁcients from a Chebyshev ﬁt polynomial coefﬁcients of a shifted polynomial inverse of chebpc; use to economize power series Pad´ approximant from power series coefﬁcients e rational ﬁt by least-squares method logarithm of gamma function factorial function binomial coefﬁcients function logarithm of factorial function

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 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).

Computer Programs by Chapter and Section

xxi

6.1 6.2 6.2 6.2 6.2 6.2 6.2 6.2 6.3 6.3 6.4 6.4 6.5 6.5 6.5 6.5 6.5 6.5 6.6 6.6 6.6 6.6 6.6 6.6 6.7 6.7 6.7 6.7 6.7 6.8 6.9 6.9 6.10 6.11 6.11 6.11 6.11 6.11 6.11 6.11 6.11 6.12 6.12 6.12 7.1 7.1

beta gammp gammq gser gcf erff erffc erfcc expint ei betai betacf bessj0 bessy0 bessj1 bessy1 bessy bessj bessi0 bessk0 bessi1 bessk1 bessk bessi bessjy beschb bessik airy sphbes plgndr frenel cisi dawson rf rd rj rc ellf elle ellpi sncndn hypgeo hypser hypdrv ran0 ran1

beta function incomplete gamma function complement of incomplete gamma function series used by gammp and gammq continued fraction used by gammp and gammq error function complementary error function complementary error function, concise routine exponential integral E n exponential integral Ei incomplete beta function continued fraction used by betai Bessel function J 0 Bessel function Y 0 Bessel function J 1 Bessel function Y 1 Bessel function Y of general integer order Bessel function J of general integer order modiﬁed Bessel function I 0 modiﬁed Bessel function K 0 modiﬁed Bessel function I 1 modiﬁed Bessel function K 1 modiﬁed Bessel function K of integer order modiﬁed Bessel function I of integer order Bessel functions of fractional order Chebyshev expansion used by bessjy modiﬁed Bessel functions of fractional order Airy functions spherical Bessel functions j n and yn Legendre polynomials, associated (spherical harmonics) Fresnel integrals S(x) and C(x) cosine and sine integrals Ci and Si Dawson’s integral Carlson’s elliptic integral of the ﬁrst kind Carlson’s elliptic integral of the second kind Carlson’s elliptic integral of the third kind Carlson’s degenerate elliptic integral Legendre elliptic integral of the ﬁrst kind Legendre elliptic integral of the second kind Legendre elliptic integral of the third kind Jacobian elliptic functions complex hypergeometric function complex hypergeometric function, series evaluation complex hypergeometric function, derivative of random deviate by Park and Miller minimal standard random deviate, minimal standard plus shufﬂe

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 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).

xxii 7.1 7.1 7.2 7.2 7.3 7.3 7.3 7.4 7.4 7.5 7.5 7.7 7.8 7.8 7.8 7.8 8.1 8.1 8.1 8.2 8.2 8.3 8.4 8.4 8.4 8.5 8.5 8.5 8.6 8.6 9.0 9.1 9.1 9.1 9.2 9.2 9.2 9.3 9.4 9.4 9.5 9.5 9.5 9.5 ran2 ran3 expdev gasdev gamdev poidev bnldev irbit1 irbit2 psdes ran4 sobseq vegas rebin miser ranpt piksrt piksr2 shell sort sort2 hpsort indexx sort3 rank select selip hpsel eclass eclazz scrsho zbrac zbrak rtbis rtflsp rtsec zriddr zbrent rtnewt rtsafe laguer zroots zrhqr qroot

Computer Programs by Chapter and Section

random deviate by L’Ecuyer long period plus shufﬂe random deviate by Knuth subtractive method exponential random deviates normally distributed random deviates gamma-law distribution random deviates Poisson distributed random deviates binomial distributed random deviates random bit sequence random bit sequence “pseudo-DES” hashing of 64 bits random deviates from DES-like hashing Sobol’s quasi-random sequence adaptive multidimensional Monte Carlo integration sample rebinning used by vegas recursive multidimensional Monte Carlo integration get random point, used by miser sort an array by straight insertion sort two arrays by straight insertion sort an array by Shell’s method sort an array by quicksort method sort two arrays by quicksort method sort an array by heapsort method construct an index for an array sort, use an index to sort 3 or more arrays construct a rank table for an array ﬁnd the N th largest in an array ﬁnd the N th largest, without altering an array ﬁnd M largest values, without altering an array determine equivalence classes from list determine equivalence classes from procedure graph a function to search for roots outward search for brackets on roots inward search for brackets on roots ﬁnd root of a function by bisection ﬁnd root of a function by false-position ﬁnd root of a function by secant method ﬁnd root of a function by Ridders’ method ﬁnd root of a function by Brent’s method ﬁnd root of a function by Newton-Raphson ﬁnd root of a function by Newton-Raphson and bisection ﬁnd a root of a polynomial by Laguerre’s method roots of a polynomial by Laguerre’s method with deﬂation roots of a polynomial by eigenvalue methods complex or double root of a polynomial, Bairstow

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 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).

Computer Programs by Chapter and Section

xxiii

9.6 9.7 9.7 9.7 9.7 9.7 10.1 10.1 10.2 10.3 10.4 10.4 10.5 10.5 10.5 10.6 10.6 10.6 10.7 10.8 10.8 10.8 10.8 10.9 10.9 10.9 10.9 10.9 10.9 10.9 10.9 11.1 11.1 11.2 11.3 11.5 11.5 11.6 12.2 12.3 12.3 12.3 12.3 12.3

mnewt lnsrch newt fdjac fmin broydn mnbrak golden brent dbrent amoeba amotry powell linmin f1dim frprmn dlinmin df1dim dfpmin simplx simp1 simp2 simp3 anneal revcst reverse trncst trnspt metrop amebsa amotsa jacobi eigsrt tred2 tqli balanc elmhes hqr four1 twofft realft sinft cosft1 cosft2

Newton’s method for systems of equations search along a line, used by newt globally convergent multi-dimensional Newton’s method ﬁnite-difference Jacobian, used by newt norm of a vector function, used by newt secant method for systems of equations bracket the minimum of a function ﬁnd minimum of a function by golden section search ﬁnd minimum of a function by Brent’s method ﬁnd minimum of a function using derivative information minimize in N -dimensions by downhill simplex method evaluate a trial point, used by amoeba minimize in N -dimensions by Powell’s method minimum of a function along a ray in N -dimensions function used by linmin minimize in N -dimensions by conjugate gradient minimum of a function along a ray using derivatives function used by dlinmin minimize in N -dimensions by variable metric method linear programming maximization of a linear function linear programming, used by simplx linear programming, used by simplx linear programming, used by simplx traveling salesman problem by simulated annealing cost of a reversal, used by anneal do a reversal, used by anneal cost of a transposition, used by anneal do a transposition, used by anneal Metropolis algorithm, used by anneal simulated annealing in continuous spaces evaluate a trial point, used by amebsa eigenvalues and eigenvectors of a symmetric matrix eigenvectors, sorts into order by eigenvalue Householder reduction of a real, symmetric matrix eigensolution of a symmetric tridiagonal matrix balance a nonsymmetric matrix reduce a general matrix to Hessenberg form eigenvalues of a Hessenberg matrix fast Fourier transform (FFT) in one dimension fast Fourier transform of two real functions fast Fourier transform of a single real function fast sine transform fast cosine transform with endpoints “staggered” fast cosine transform

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 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).

xxiv 12.4 12.5 12.6 12.6 13.1 13.2 13.4 13.6 13.6 13.6 13.7 13.8 13.8 13.8 13.9 13.9 13.10 13.10 13.10 13.10 13.10 14.1 14.2 14.2 14.2 14.2 14.2 14.3 14.3 14.3 14.3 14.3 14.4 14.4 14.5 14.6 14.6 14.6 14.6 14.7 14.7 14.7 14.7 14.8 fourn rlft3 fourfs fourew convlv correl spctrm memcof fixrts predic evlmem period fasper spread dftcor dftint wt1 daub4 pwtset pwt wtn moment ttest avevar tutest tptest ftest chsone chstwo ksone kstwo probks cntab1 cntab2 pearsn spear crank kendl1 kendl2 ks2d1s quadct quadvl ks2d2s savgol

Computer Programs by Chapter and Section

fast Fourier transform in multidimensions FFT of real data in two or three dimensions FFT for huge data sets on external media rewind and permute ﬁles, used by fourfs convolution or deconvolution of data using FFT correlation or autocorrelation of data using FFT power spectrum estimation using FFT evaluate maximum entropy (MEM) coefﬁcients reﬂect roots of a polynomial into unit circle linear prediction using MEM coefﬁcients power spectral estimation from MEM coefﬁcients power spectrum of unevenly sampled data power spectrum of unevenly sampled larger data sets extirpolate value into array, used by fasper compute endpoint corrections for Fourier integrals high-accuracy Fourier integrals one-dimensional discrete wavelet transform Daubechies 4-coefﬁcient wavelet ﬁlter initialize coefﬁcients for pwt partial wavelet transform multidimensional discrete wavelet transform calculate moments of a data set Student’s t-test for difference of means calculate mean and variance of a data set Student’s t-test for means, case of unequal variances Student’s t-test for means, case of paired data F -test for difference of variances chi-square test for difference between data and model chi-square test for difference between two data sets Kolmogorov-Smirnov test of data against model Kolmogorov-Smirnov test between two data sets Kolmogorov-Smirnov probability function contingency table analysis using chi-square contingency table analysis using entropy measure Pearson’s correlation between two data sets Spearman’s rank correlation between two data sets replaces array elements by their rank correlation between two data sets, Kendall’s tau contingency table analysis using Kendall’s tau K–S test in two dimensions, data vs. model count points by quadrants, used by ks2d1s quadrant probabilities, used by ks2d1s K–S test in two dimensions, data vs. data Savitzky-Golay smoothing coefﬁcients
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 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).

Computer Programs by Chapter and Section

xxv

15.2 15.3 15.3 15.4 15.4 15.4 15.4 15.4 15.4 15.5 15.5 15.5 15.7 15.7 16.1 16.1 16.2 16.2 16.2 16.3 16.4 16.4 16.4 16.5 16.6 16.6 16.6 16.6 16.6 17.1 17.2 17.3 17.3 17.3 17.3 17.4 17.4 17.4 17.4 18.1 18.1 18.2 18.3 18.3

fit fitexy chixy lfit covsrt svdfit svdvar fpoly fleg mrqmin mrqcof fgauss medfit rofunc rk4 rkdumb rkqs rkck odeint mmid bsstep pzextr rzextr stoerm stiff jacobn derivs simpr stifbs shoot shootf solvde bksub pinvs red sfroid difeq sphoot sphfpt fred2 fredin voltra wwghts kermom

least-squares ﬁt data to a straight line ﬁt data to a straight line, errors in both x and y used by fitexy to calculate a χ 2 general linear least-squares ﬁt by normal equations rearrange covariance matrix, used by lfit linear least-squares ﬁt by singular value decomposition variances from singular value decomposition ﬁt a polynomial using lfit or svdfit ﬁt a Legendre polynomial using lfit or svdfit nonlinear least-squares ﬁt, Marquardt’s method used by mrqmin to evaluate coefﬁcients ﬁt a sum of Gaussians using mrqmin ﬁt data to a straight line robustly, least absolute deviation ﬁt data robustly, used by medfit integrate one step of ODEs, fourth-order Runge-Kutta integrate ODEs by fourth-order Runge-Kutta integrate one step of ODEs with accuracy monitoring Cash-Karp-Runge-Kutta step used by rkqs integrate ODEs with accuracy monitoring integrate ODEs by modiﬁed midpoint method integrate ODEs, Bulirsch-Stoer step polynomial extrapolation, used by bsstep rational function extrapolation, used by bsstep integrate conservative second-order ODEs integrate stiff ODEs by fourth-order Rosenbrock sample Jacobian routine for stiff sample derivatives routine for stiff integrate stiff ODEs by semi-implicit midpoint rule integrate stiff ODEs, Bulirsch-Stoer step solve two point boundary value problem by shooting ditto, by shooting to a ﬁtting point two point boundary value problem, solve by relaxation backsubstitution, used by solvde diagonalize a sub-block, used by solvde reduce columns of a matrix, used by solvde spheroidal functions by method of solvde spheroidal matrix coefﬁcients, used by sfroid spheroidal functions by method of shoot spheroidal functions by method of shootf solve linear Fredholm equations of the second kind interpolate solutions obtained with fred2 linear Volterra equations of the second kind quadrature weights for an arbitrarily singular kernel sample routine for moments of a singular kernel

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 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).

xxvi 18.3 18.3 19.5 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 19.6 20.1 20.2 20.3 20.3 20.3 20.4 20.4 20.4 20.4 20.5 20.5 20.5 20.6 20.6 20.6 20.6 20.6 20.6 20.6 quadmx fredex sor mglin rstrct interp addint slvsml relax resid copy fill0 mgfas relax2 slvsm2 lop matadd matsub anorm2 machar igray icrc1 icrc decchk hufmak hufapp hufenc hufdec arcmak arcode arcsum mpops mpmul mpinv mpdiv mpsqrt mp2dfr mppi

Computer Programs by Chapter and Section

sample routine for a quadrature matrix example of solving a singular Fredholm equation elliptic PDE solved by successive overrelaxation method linear elliptic PDE solved by multigrid method half-weighting restriction, used by mglin, mgfas bilinear prolongation, used by mglin, mgfas interpolate and add, used by mglin solve on coarsest grid, used by mglin Gauss-Seidel relaxation, used by mglin calculate residual, used by mglin utility used by mglin, mgfas utility used by mglin nonlinear elliptic PDE solved by multigrid method Gauss-Seidel relaxation, used by mgfas solve on coarsest grid, used by mgfas applies nonlinear operator, used by mgfas utility used by mgfas utility used by mgfas utility used by mgfas diagnose computer’s ﬂoating arithmetic Gray code and its inverse cyclic redundancy checksum, used by icrc cyclic redundancy checksum decimal check digit calculation or veriﬁcation construct a Huffman code append bits to a Huffman code, used by hufmak use Huffman code to encode and compress a character use Huffman code to decode and decompress a character construct an arithmetic code encode or decode a character using arithmetic coding add integer to byte string, used by arcode multiple precision arithmetic, simpler operations multiple precision multiply, using FFT methods multiple precision reciprocal multiple precision divide and remainder multiple precision square root multiple precision conversion to decimal base multiple precision example, compute many digits of π

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 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).

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 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).

Chapter 1.
1.0 Introduction

Preliminaries

This book, like its predecessor edition, is supposed to teach you methods of numerical computing that are practical, efﬁcient, and (insofar as possible) elegant. We presume throughout this book that you, the reader, have particular tasks that you want to get done. We view our job as educating you on how to proceed. Occasionally we may try to reroute you brieﬂy onto a particularly beautiful side road; but by and large, we will guide you along main highways that lead to practical destinations. Throughout this book, you will ﬁnd us fearlessly editorializing, telling you what you should and shouldn’t do. This prescriptive tone results from a conscious decision on our part, and we hope that you will not ﬁnd it irritating. We do not claim that our advice is infallible! Rather, we are reacting against a tendency, in the textbook literature of computation, to discuss every possible method that has ever been invented, without ever offering a practical judgment on relative merit. We do, therefore, offer you our practical judgments whenever we can. As you gain experience, 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 being the language of this version of Numerical Recipes (Second Edition). The book Numerical Recipes in FORTRAN (Second Edition) is separately available, if you prefer to program in that language. Earlier editions of Numerical Recipes in Pascal and Numerical Recipes Routines and Examples in BASIC are also available; while not containing the additional material of the Second Edition versions in C and FORTRAN, these versions are perfectly serviceable if Pascal or BASIC is your language of choice. When we include programs in the text, they look like this:
#include #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 explaining their calling sequence. This routine calculates the phases of the moon. Given an integer n and a code nph for the phase desired (nph = 0 for new moon, 1 for ﬁrst quarter, 2 for full, 3 for last quarter), the routine returns the Julian Day Number jd, and the fractional part of a day frac to 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; c=n+nph/4.0; This is how we comment an individual line.

1

2

Chapter 1.

Preliminaries

}

t=c/1236.85; t2=t*t; as=359.2242+29.105356*c; You aren’t really intended to understand am=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 error i=(int)(xtra >= 0.0 ? floor(xtra) : ceil(xtra-1.0)); conditions. *jd += i; *frac=xtra-i;

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 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 the syntax of the function deﬁnition above looks strange to you, then you are probably used to the older Kernighan and Ritchie (“K&R”) syntax, rather than that of the newer ANSI C. In this edition, we adopt ANSI C as our standard. You might want to 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 statement like nrerror("some error message");. The function nrerror() is part of a small ﬁle of utility programs, nrutil.c, listed in Appendix B at the back of the book. This Appendix includes a number of other utilities that we will describe later in this chapter. Function nrerror() prints the indicated error message to your stderr device (usually your terminal screen), and then invokes the function exit(), which terminates execution. The function exit() is in every C library we know of; but if you ﬁnd it missing, you can modify nrerror() so that it does anything else that will halt execution. For example, you can have it pause for input from the keyboard, and then manually interrupt execution. In some applications, you will want to modify nrerror() to do more sophisticated error handling, for example to transfer control somewhere else, with an error ﬂag or error code set. We will have more to say about the C programming language, its conventions and style, in §1.1 and §1.2.

Computational Environment and Program Validation
Our goal is that the programs in this book be as portable as possible, across different platforms (models of computer), across different operating systems, and across different C compilers. C was designed with this type of portability in mind. Nevertheless, we have found that there is no substitute for actually checking all programs on a variety of compilers, in the process uncovering differences in library structure or contents, and even occasional differences in allowed syntax. As surrogates for the large number of possible combinations, we have tested all the programs in this book on the combinations of machines, operating systems, and compilers shown on the accompanying table. More generally, the programs should run without modiﬁcation on any compiler that implements the ANSI C standard, as described for example in Harbison and Steele’s excellent book [1]. With small modiﬁcations, 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 to watch out for is that ANSI C requires the memory allocation functions malloc()

1.0 Introduction

3

Tested Machines and Compilers Hardware IBM PC compatible 486/33 IBM PC compatible 486/33 IBM RS/6000 DECstation 5000/25 DECsystem 5400 Sun SPARCstation 2 DECstation 5000/200 Sun SPARCstation 2 O/S Version MS-DOS 5.0/Windows 3.1 MS-DOS 5.0 AIX 3.2 ULTRIX 4.2A ULTRIX 4.1 SunOS 4.1 ULTRIX 4.2 SunOS 4.1 Compiler Version Microsoft C/C++ 7.0 Borland C/C++ 2.0 IBM xlc 1.02 CodeCenter (Saber) C 3.1.1 GNU C Compiler 2.1 GNU C Compiler 1.40 DEC RISC C 2.1* Sun cc 1.1*
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 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).

*compiler version does not fully implement ANSI C; only K&R validated

and free() to be declared via the header stdlib.h; some older compilers require them to be declared with the header ﬁle malloc.h, while others regard them as inherent in the language and require no header ﬁle at all. In validating the programs, we have taken the program source code directly from the machine-readable form of the book’s manuscript, to decrease the chance of propagating typographical errors. “Driver” or demonstration programs that we used as part of our validations are available separately as the Numerical Recipes Example Book (C), as well as in machine-readable form. If you plan to use more than a few of the programs in this book, or if you plan to use programs in this book on more than one different computer, then you may ﬁnd it useful to obtain a copy of 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 beneﬁtted from the experience of the many readers who have written to us. If you ﬁnd a new bug, please document it and tell us!

Compatibility with the First Edition
If you are accustomed to the Numerical Recipes routines of the First Edition, rest assured: 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 you will soon become equally familiar with the added capabilities of the more than 100 routines that are new to this edition. We have retired a small number of First Edition routines, those that we believe to 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 to both editions have alterations in their calling interfaces, so are not directly “plug compatible.” 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 taken for the comparison). If you have written software of any appreciable complexity

4

Chapter 1.

Preliminaries

Previous Routines Omitted from This Edition Name(s) adi cosft cel, el2 des, desks mdian1, mdian2 qcksrt rkqc smooft sparse Replacement(s) mglin or mgfas cosft1 or cosft2 rf, rd, rj, rc ran4 now uses psdes select, selip sort rkqs Comment better method choice of boundary conditions better algorithms was too slow more general name change (sort is now hpsort) better method
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 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).

use convlv with coefﬁcients from savgol linbcg more general

that is dependent on First Edition routines, we do not recommend blindly replacing them by the corresponding routines in this book. We do recommend that any new programming efforts use the new routines.

About References
You will ﬁnd references, and suggestions for further reading, listed at the end of most sections of this book. References are cited in the text by bracketed numbers like this [3]. Because computer algorithms often circulate informally for quite some time before appearing in a published form, the task of uncovering “primary literature” is sometimes quite difﬁcult. We have not attempted this, and we do not pretend to any degree of bibliographical completeness in this book. For topics where a substantial secondary literature exists (discussion in textbooks, reviews, etc.) we have consciously limited our references to a few of the more useful secondary sources, especially those with good references to the primary literature. Where the existing secondary literature is insufﬁcient, we give references to a few primary sources that are intended to serve as starting points for further reading, not as complete bibliographies for the ﬁeld. The order in which references are listed is not necessarily signiﬁcant. It reﬂects a compromise between listing cited references in the order cited, and listing suggestions for further reading in a roughly prioritized order, with the most useful ones ﬁrst. The remaining three sections of this chapter review some basic concepts of programming (control structures, etc.), discuss a set of conventions speciﬁc to C that we have adopted in this book, and introduce some fundamental concepts in numerical analysis (roundoff error, etc.). Thereafter, we plunge into the substantive material of the book.
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). [1]

1.1 Program Organization and Control Structures

5

Kernighan, B., and Ritchie, D. 1978, The C Programming Language (Englewood Cliffs, NJ: Prentice-Hall). [2] [Reference for K&R “traditional” C. Later editions of this book conform to the ANSI C standard.] Meeus, J. 1982, Astronomical Formulae for Calculators, 2nd ed., revised and enlarged (Richmond, VA: Willmann-Bell). [3]
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 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).

1.1 Program Organization and Control Structures
We sometimes like to point out the close analogies between computer programs, on the one hand, and written poetry or written musical scores, on the other. All three present themselves as visual media, symbols on a two-dimensional page or computer screen. Yet, in all three cases, the visual, two-dimensional, frozen-in-time representation communicates (or is supposed to communicate) something rather different, namely a process that unfolds in time. A poem is meant to be read; music, played; a program, executed as a sequential series of computer instructions. In all three cases, the target of the communication, in its visual form, is a human being. The goal is to transfer to him/her, as efﬁciently as can be accomplished, the greatest degree of understanding, in advance, of how the process will unfold in time. In poetry, this human target is the reader. In music, it is the performer. In programming, it is the program user. Now, you may object that the target of communication of a program is not a human but a computer, that the program user is only an irrelevant intermediary, a lackey who feeds the machine. This is perhaps the case in the situation where the business executive pops a diskette into a desktop computer and feeds that computer a black-box program in binary executable form. The computer, in this case, doesn’t much care whether that program was written with “good programming practice” or not. We envision, however, that you, the readers of this book, are in quite a different situation. You need, or want, to know not just what a program does, but also how it does it, so that you can tinker with it and modify it to your particular application. You need others to be able to see what you have done, so that they can criticize or admire. In such cases, where the desired goal is maintainable or reusable code, the targets of a program’s communication are surely human, not machine. One key to achieving good programming practice is to recognize that programming, music, and poetry — all three being symbolic constructs of the human brain — are naturally structured into hierarchies that have many different nested levels. Sounds (phonemes) form small meaningful units (morphemes) which in turn form words; words group into phrases, which group into sentences; sentences make paragraphs, and these are organized into higher levels of meaning. Notes form musical phrases, which form themes, counterpoints, harmonies, etc.; which form movements, which form concertos, symphonies, and so on. The structure in programs is equally hierarchical. Appropriately, good programming practice brings different techniques to bear on the different levels [1-3]. At a low level is the ascii character set. Then, constants, identiﬁers, operands, operators.

6

Chapter 1.

Preliminaries

Then program statements, like a[j+1]=b+c/3.0;. Here, the best programming advice is simply be clear, or (correspondingly) don’t be too tricky. You might momentarily be proud of yourself at writing the single line k=(2-j)*(1+3*j)/2; if you want to permute cyclically one of the values j = (0, 1, 2) into respectively k = (1, 2, 0). You will regret it later, however, when you try to understand that line. Better, and likely also faster, is k=j+1; if (k == 3) k=0;

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 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).

Many programming stylists would even argue for the ploddingly literal switch (j) { case 0: k=1; break; case 1: k=2; break; case 2: k=0; break; default: { fprintf(stderr,"unexpected value for j"); exit(1); } }

on the grounds that it is both clear and additionally safeguarded from wrong assumptions about the possible values of j. Our preference among the implementations is for the middle one. In this simple example, we have in fact traversed several levels of hierarchy: Statements frequently come in “groups” or “blocks” which make sense only taken as a whole. The middle fragment above is one example. Another is swap=a[j]; a[j]=b[j]; b[j]=swap;

which makes immediate sense to any programmer as the exchange of two variables, while ans=sum=0.0; n=1;

is very likely to be an initialization of variables prior to some iterative process. This level of hierarchy in a program is usually evident to the eye. It is good programming practice to put in comments at this level, e.g., “initialize” or “exchange variables.” The next level is that of control structures. These are things like the switch construction in the example above, for loops, and so on. This level is sufﬁciently important, and relevant to the hierarchical level of the routines in this book, that we will come back to it just below. At still higher levels in the hierarchy, we have functions and modules, and the whole “global” organization of the computational task to be done. In the musical analogy, we are now at the level of movements and complete works. At these levels,

1.1 Program Organization and Control Structures

7

modularization and encapsulation become important programming concepts, the general idea being that program units should interact with one another only through clearly deﬁned and narrowly circumscribed interfaces. Good modularization practice is an essential prerequisite to the success of large, complicated software projects, especially those employing the efforts of more than one programmer. It is also good practice (if not quite as essential) in the less massive programming tasks that an individual scientist, or reader of this book, encounters. Some computer languages, such as Modula-2 and C++, promote good modularization with higher-level language constructs absent in C. In Modula-2, for example, functions, type deﬁnitions, and data structures can be encapsulated into “modules” that communicate through declared public interfaces and whose internal workings are hidden from the rest of the program [4]. In the C++ language, the key concept is “class,” a user-deﬁnable generalization of data type that provides for data hiding, automatic initialization of data, memory management, dynamic typing, and operator overloading (i.e., the user-deﬁnable extension of operators like + and * so as to be appropriate to operands in any particular class) [5]. Properly used in deﬁning the data structures that are passed between program units, classes can clarify and circumscribe these units’ public interfaces, reducing the chances of programming error and also allowing a considerable degree of compile-time and run-time error checking. Beyond modularization, though depending on it, lie the concepts of objectoriented programming. Here a programming language, such as C++ or Turbo Pascal 5.5 [6], allows a module’s public interface to accept redeﬁnitions of types or actions, and these redeﬁnitions become shared all the way down through the module’s hierarchy (so-called polymorphism). For example, a routine written to invert a matrix of real numbers could — dynamically, at run time — be made able to handle complex numbers by overloading complex data types and corresponding deﬁnitions of the arithmetic operations. Additional concepts of inheritance (the ability to deﬁne a data type that “inherits” all the structure of another type, plus additional structure of its own), and object extensibility (the ability to add functionality to a module without access to its source code, e.g., at run time), also come into play. We have not attempted to modularize, or make objects out of, the routines in this book, for at least two reasons. First, the chosen language, C, does not really make this possible. Second, we envision that you, the reader, might want to incorporate the algorithms in this book, a few at a time, into modules or objects with a structure of your own choosing. There does not exist, at present, a standard or accepted set of “classes” for scientiﬁc object-oriented computing. While we might have tried to invent such a set, doing so would have inevitably tied the algorithmic content of the ˆ book (which is its raison d’ etre) to some rather speciﬁc, and perhaps haphazard, set of choices regarding class deﬁnitions. On the other hand, we are not unfriendly to the goals of modular and objectoriented programming. Within the limits of C, we have therefore tried to structure our programs to be “object friendly.” That is one reason we have adopted ANSI C with its function prototyping as our default C dialect (see §1.2). Also, within our implementation sections, we have paid particular attention to the practices of structured programming, as we now discuss.

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 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).

8

Chapter 1.

Preliminaries

Control Structures
An executing program unfolds in time, but not strictly in the linear order in which the statements are written. Program statements that affect the order in which statements are executed, or that affect whether statements are executed, are called control statements. Control statements never make useful sense by themselves. They make sense only in the context of the groups or blocks of statements that they in turn control. If you think of those blocks as paragraphs containing sentences, then the control statements are perhaps best thought of as the indentation of the paragraph and the punctuation between the sentences, not the words within the sentences. We can now say what the goal of structured programming is. It is to make program control manifestly apparent in the visual presentation of the program. You see that this goal has nothing at all to do with how the computer sees the program. As already remarked, computers don’t care whether you use structured programming or not. Human readers, however, do care. You yourself will also care, once you discover how much easier it is to perfect and debug a well-structured program than one whose control structure is obscure. You accomplish the goals of structured programming in two complementary ways. First, you acquaint yourself with the small number of essential control structures that occur over and over again in programming, and that are therefore given convenient representations in most programming languages. You should learn to think about your programming tasks, insofar as possible, exclusively in terms of these standard control structures. In writing programs, you should get into the habit of representing these standard control structures in consistent, conventional ways. “Doesn’t this inhibit creativity?” our students sometimes ask. Yes, just as Mozart’s creativity was inhibited by the sonata form, or Shakespeare’s by the metrical requirements of the sonnet. The point is that creativity, when it is meant to communicate, does well under the inhibitions of appropriate restrictions on format. Second, you avoid, insofar as possible, control statements whose controlled blocks or objects are difﬁcult to discern at a glance. This means, in practice, that you must try to avoid named labels on statements and goto’s. It is not the goto’s that are dangerous (although they do interrupt one’s reading of a program); the named statement labels are the hazard. In fact, whenever you encounter a named statement label while reading a program, you will soon become conditioned to get a sinking feeling in the pit of your stomach. Why? Because the following questions will, by habit, immediately spring to mind: Where did control come from in a branch to this label? It could be anywhere in the routine! What circumstances resulted in a branch to this label? They could be anything! Certainty becomes uncertainty, understanding dissolves into a morass of possibilities. Some examples are now in order to make these considerations more concrete (see Figure 1.1.1).

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 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).

Catalog of Standard Structures
Iteration. In C, simple iteration is performed with a for loop, for example

for (j=2;j 3) if (a > 3) b += 1; else b -= 1;

/* questionable! */

As judged by the indentation used on successive lines, the intent of the writer of this code is the following: ‘If b is greater than 3 and a is greater than 3, then increment b. If b is not greater than 3, then decrement b.’ According to the rules of C, however, the actual meaning is ‘If b is greater than 3, then evaluate a. If a is greater than 3, then increment b, and if a is less than or equal to 3, decrement b.’ The point is that an else clause is associated with the most recent open if statement, no matter how you lay it out on the page. Such confusions in meaning are easily resolved by the inclusion of braces. They may in some instances be technically superﬂuous; nevertheless, they clarify your intent and improve the program. The above fragment should be written as if (b > 3) { if (a > 3) b += 1; } else { b -= 1; }

Here is a working program that consists dominantly of if control statements:
#include #define IGREG (15+31L*(10+12L*1582)) Gregorian Calendar adopted Oct. 15, 1582.

long julday(int mm, int id, int iyyy) In this routine julday returns the Julian Day Number that begins at noon of the calendar date speciﬁed by month mm, day id, and year iyyy, all integer variables. Positive year signiﬁes A.D.; negative, B.C. Remember that the year after 1 B.C. was 1 A.D. { void nrerror(char error_text[]);

12 long jul; int ja,jy=iyyy,jm;

Chapter 1.

Preliminaries

}

if (jy == 0) nrerror("julday: there is no year zero."); if (jy < 0) ++jy; if (mm > 2) { Here is an example of a block IF-structure. jm=mm+1; } else { --jy; jm=mm+13; } jul = (long) (floor(365.25*jy)+floor(30.6001*jm)+id+1720995); if (id+31L*(mm+12L*iyyy) >= IGREG) { Test whether to change to Gregorian Calja=(int)(0.01*jy); endar. jul += 2-ja+(int) (0.25*ja); } return jul;

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 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).

(Astronomers number each 24-hour period, starting and ending at noon, with a unique integer, the Julian Day Number [7]. Julian Day Zero was a very long time ago; a convenient reference point is that Julian Day 2440000 began at noon of May 23, 1968. If you know the Julian Day Number that begins at noon of a given calendar date, then the day of the week of that date is obtained by adding 1 and taking the result modulo base 7; a zero answer corresponds to Sunday, 1 to Monday, . . . , 6 to Saturday.) While iteration. Most languages (though not FORTRAN, incidentally) provide for structures like the following C example: while (n < 1000) { n *= 2; j += 1; }

It is the particular feature of this structure that the control-clause (in this case n < 1000) is evaluated before each iteration. If the clause is not true, the enclosed statements will not be executed. In particular, if this code is encountered at a time when n is greater than or equal to 1000, the statements will not even be executed once. Do-While iteration. Companion to the while iteration is a related controlstructure that tests its control-clause at the end of each iteration. In C, it looks like this: do { n *= 2; j += 1; } while (n < 1000);

In this case, the enclosed statements will be executed at least once, independent of the initial value of n. Break. In this case, you have a loop that is repeated indeﬁnitely until some condition tested somewhere in the middle of the loop (and possibly tested in more

1.1 Program Organization and Control Structures

13

than one place) becomes true. At that point you wish to exit the loop and proceed with what comes after it. In C the structure is implemented with the simple break statement, which terminates execution of the innermost for, while, do, or switch construction and proceeds to the next sequential instruction. (In Pascal and standard FORTRAN, this structure requires the use of statement labels, to the detriment of clear programming.) A typical usage of the break statement is: for(;;) { [statements before the test] if (...) break; [statements after the test] } [next sequential instruction]

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 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).

Here is a program that uses several different iteration structures. One of us was once asked, for a scavenger hunt, to ﬁnd the date of a Friday the 13th on which the moon was full. This is a program which accomplishes that task, giving incidentally all other Fridays the 13th as a by-product.
#include #include #define ZON -5.0 #define IYBEG 1900 #define IYEND 2000

Time zone −5 is Eastern Standard Time. The range of dates to be searched.

int main(void) /* Program badluk */ { void flmoon(int n, int nph, long *jd, float *frac); long julday(int mm, int id, int iyyy); int ic,icon,idwk,im,iyyy,n; float timzon = ZON/24.0,frac; long jd,jday; printf("\nFull moons on Friday the 13th from %5d to %5d\n",IYBEG,IYEND); for (iyyy=IYBEG;iyyy= jd ? 1 : -1); if (ic == (-icon)) break; Another break, case of no match. icon=ic; n += ic; }

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 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 you are merely curious, there were (or will be) occurrences of a full moon on Friday the 13th (time zone GMT−5) on: 3/13/1903, 10/13/1905, 6/13/1919, 1/13/1922, 11/13/1970, 2/13/1987, 10/13/2000, 9/13/2019, and 8/13/2049. Our advice is to avoid them. Every proOther “standard” structures. gramming language has some number of “goodies” that the designer just couldn’t resist throwing in. They seemed like a good idea at the time. Unfortunately they don’t stand the test of time! Your program becomes difﬁcult to translate into other languages, and difﬁcult to read (because rarely used structures are unfamiliar to the reader). You can almost always accomplish the supposed conveniences of these structures in other ways. In C, the most problematic control structure is the switch...case...default construction (see Figure 1.1.1), which has historically been burdened by uncertainty, from compiler to compiler, about what data types are allowed in its control expression. Data types char and int are universally supported. For other data types, e.g., float or double, the structure should be replaced by a more recognizable and translatable if. . .else construction. ANSI C allows the control expression to be of type long, but many older compilers do not. The continue; construction, while benign, can generally be replaced by an if construction with no loss of clarity.

About “Advanced Topics”
Material set in smaller type, like this, signals an “advanced topic,” either one outside of the main argument of the chapter, or else one requiring of you more than the usual assumed mathematical background, or else (in a few cases) a discussion that is more speculative or an algorithm that is less well-tested. Nothing important will be lost if you skip the advanced topics on a ﬁrst reading of the book. You may have noticed that, by its looping over the months and years, the program badluk avoids using any algorithm for converting a Julian Day Number back into a calendar date. A routine for doing just this is not very interesting structurally, but it is occasionally useful:
#include #define IGREG 2299161 void caldat(long julian, int *mm, int *id, int *iyyy) Inverse of the function julday given above. Here julian is input as a Julian Day Number, and the routine outputs mm,id, and iyyy as the month, day, and year on which the speciﬁed Julian Day started at noon. { long ja,jalpha,jb,jc,jd,je;

1.2 Some C Conventions for Scientiﬁc Computing

15

}

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 < 0) { Make day number positive by adding integer number of ja=julian+36525*(1-julian/36525); Julian centuries, then subtract them oﬀ } else at the end. ja=julian; jb=ja+1524; jc=(long)(6680.0+((double) (jb-2439870)-122.1)/365.25); jd=(long)(365*jc+(0.25*jc)); je=(long)((jb-jd)/30.6001); *id=jb-jd-(long) (30.6001*je); *mm=je-1; if (*mm > 12) *mm -= 12; *iyyy=jc-4715; if (*mm > 2) --(*iyyy); if (*iyyy < > = == != & ^ | && || ? : function call array element structure or union member pointer reference to structure logical not bitwise complement unary minus increment decrement address of contents of cast to type size in bytes multiply divide remainder add subtract bitwise left shift bitwise right shift arithmetic arithmetic arithmetic arithmetic less than greater than less than or equal to greater than or equal to left-to-right

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 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).

right-to-left

left-to-right

left-to-right left-to-right left-to-right

arithmetic equal arithmetic not equal bitwise and bitwise exclusive or bitwise or logical and logical or conditional expression

left-to-right left-to-right left-to-right left-to-right left-to-right left-to-right right-to-left right-to-left

= assignment operator also += -= *= /= %= = &= ^= |= , sequential expression

left-to-right

We have already alluded to the problem of computing small integer powers of numbers, most notably the square and cube. The omission of this operation from C is perhaps the language’s most galling insult to the scientiﬁc programmer. All good FORTRAN compilers recognize expressions like (A+B)**4 and produce in-line code, in this case with only one add and two multiplies. It is typical for constant integer powers up to 12 to be thus recognized.

1.2 Some C Conventions for Scientiﬁc Computing

27

In C, the mere problem of squaring is hard enough! Some people “macro-ize” the operation as
#define SQR(a) ((a)*(a))

However, this is likely to produce code where SQR(sin(x)) results in two calls to the sine routine! You might be tempted to avoid this by storing the argument of the squaring function in a temporary variable: static float sqrarg; #define SQR(a) (sqrarg=(a),sqrarg*sqrarg)

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 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).

The global variable sqrarg now has (and needs to keep) scope over the whole module, which is a little dangerous. Also, one needs a completely different macro to square expressions of type int. More seriously, this macro can fail if there are two SQR operations in a single expression. Since in C the order of evaluation of pieces of the expression is at the compiler’s discretion, the value of sqrarg in one evaluation of SQR can be that from the other evaluation in the same expression, producing nonsensical results. When we need a guaranteed-correct SQR macro, we use the following, which exploits the guaranteed complete evaluation of subexpressions in a conditional expression: static float sqrarg; #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

A collection of macros for other simple operations is included in the ﬁle nrutil.h (see Appendix B) and used by many of our programs. Here are the synopses:
SQR(a) DSQR(a) FMAX(a,b) FMIN(a,b) DMAX(a,b) DMIN(a,b) IMAX(a,b) IMIN(a,b) LMAX(a,b) LMIN(a,b) SIGN(a,b) Square a float value. Square a double value. Maximum of two float values. Minimum of two float values. Maximum of two double values. Minimum of two double values. Maximum of two int values. Minimum of two int values. Maximum of two long values. Minimum of two long values. Magnitude of a times sign of b.

Scientiﬁc programming in C may someday become a bed of roses; for now, watch out for the thorns!
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). [1] AT&T Bell Laboratories 1985, The C Programmer’s Handbook (Englewood Cliffs, NJ: PrenticeHall). Kernighan, B., and Ritchie, D. 1978, The C Programming Language (Englewood Cliffs, NJ: Prentice-Hall). [Reference for K&R “traditional” C. Later editions of this book conform to the ANSI C standard.] Hogan, T. 1984, The C Programmer’s Handbook (Bowie, MD: Brady Communications).

28

Chapter 1.

Preliminaries

1.3 Error, Accuracy, and Stability
Although we assume no prior training of the reader in formal numerical analysis, we will need to presume a common understanding of a few key concepts. We will deﬁne these brieﬂy in this section. Computers store numbers not with inﬁnite precision but rather in some approximation that can be packed into a ﬁxed number of bits (binary digits) or bytes (groups of 8 bits). Almost all computers allow the programmer a choice among several different such representations or data types. Data types can differ in the number of bits utilized (the wordlength), but also in the more fundamental respect of whether the stored number is represented in ﬁxed-point (int or long) or ﬂoating-point (float or double) format. A number in integer representation is exact. Arithmetic between numbers in integer representation is also exact, with the provisos that (i) the answer is not outside the range of (usually, signed) integers that can be represented, and (ii) that division is interpreted as producing an integer result, throwing away any integer remainder. In ﬂoating-point representation, a number is represented internally by a sign bit s (interpreted as plus or minus), an exact integer exponent e, and an exact positive integer mantissa M . Taken together these represent the number s × M × B e−E (1.3.1)

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 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).

where B is the base of the representation (usually B = 2, but sometimes B = 16), and E is the bias of the exponent, a ﬁxed integer constant for any given machine and representation. An example is shown in Figure 1.3.1. Several ﬂoating-point bit patterns can represent the same number. If B = 2, for example, a mantissa with leading (high-order) zero bits can be left-shifted, i.e., multiplied by a power of 2, if the exponent is decreased by a compensating amount. Bit patterns that are “as left-shifted as they can be” are termed normalized. Most computers always produce normalized results, since these don’t waste any bits of the mantissa and thus allow a greater accuracy of the representation. Since the high-order bit of a properly normalized mantissa (when B = 2) is always one, some computers don’t store this bit at all, giving one extra bit of signiﬁcance. Arithmetic among numbers in ﬂoating-point representation is not exact, even if the operands happen to be exactly represented (i.e., have exact values in the form of equation 1.3.1). For example, two ﬂoating numbers are added by ﬁrst right-shifting (dividing by two) the mantissa of the smaller (in magnitude) one, simultaneously increasing its exponent, until the two operands have the same exponent. Low-order (least signiﬁcant) bits of the smaller operand are lost by this shifting. If the two operands differ too greatly in magnitude, then the smaller operand is effectively replaced by zero, since it is right-shifted to oblivion. The smallest (in magnitude) ﬂoating-point number which, when added to the ﬂoating-point number 1.0, produces a ﬂoating-point result different from 1.0 is termed the machine accuracy m . A typical computer with B = 2 and a 32-bit wordlength has m around 3 × 10 −8 . (A more detailed discussion of machine characteristics, and a program to determine them, is given in §20.1.) Roughly

1.3 Error, Accuracy, and Stability

29

th be is b “p it c ha ou nt ld om ”

on

xp

bi te

sig n

1⁄ 2

= 0 = 0 = 0 = 0 = 0

23

8-

-b

it m

bi t

an

tis sa

en

t

10000000 10000010 01111111 01101001 10000010 10000010

10000000000000000000000 11000000000000000000000 10000000000000000000000 1 1 1 0 ... 0 1 1 0 1 0 1 1 1 1 1 1 1 0 0 1 0 1 0 00000000000000000000000 11000000000000000000000

(a) (b) (c) (d) (e) (f )

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 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).

3 = 0
1⁄ 4

10 − 7 3+ 10 − 7

Figure 1.3.1. Floating point representations of numbers in a typical 32-bit (4-byte) format. (a) The number 1/2 (note the bias in the exponent); (b) the number 3; (c) the number 1/4; (d) the number 10−7 , represented to machine accuracy; (e) the same number 10−7 , but shifted so as to have the same −7 exponent as the number 3; with this shifting, all signiﬁcance is lost and 10 becomes zero; shifting to −7 a common exponent must occur before two numbers can be added; (f) sum of the numbers 3 + 10 , which equals 3 to machine accuracy. Even though 10−7 can be represented accurately by itself, it cannot accurately be added to a much larger number.

speaking, the machine accuracy m is the fractional accuracy to which ﬂoating-point numbers are represented, corresponding to a change of one in the least signiﬁcant bit of the mantissa. Pretty much any arithmetic operation among ﬂoating numbers should be thought of as introducing an additional fractional error of at least m . This type of error is called roundoff error. It is important to understand that m is not the smallest ﬂoating-point number that can be represented on a machine. That number depends on how many bits there are in the exponent, while m depends on how many bits there are in the mantissa. Roundoff errors accumulate with increasing amounts of calculation. If, in the course of obtaining a calculated value, you perform N such arithmetic operations, √ you might be so lucky as to have a total roundoff error on the order of N m , if the roundoff errors come in randomly up or down. (The square root comes from a random-walk.) However, this estimate can be very badly off the mark for two reasons: (i) It very frequently happens that the regularities of your calculation, or the peculiarities of your computer, cause the roundoff errors to accumulate preferentially in one direction. In this case the total will be of order N m . (ii) Some especially unfavorable occurrences can vastly increase the roundoff error of single operations. Generally these can be traced to the subtraction of two very nearly equal numbers, giving a result whose only signiﬁcant bits are those (few) low-order ones in which the operands differed. You might think that such a “coincidental” subtraction is unlikely to occur. Not always so. Some mathematical expressions magnify its probability of occurrence tremendously. For example, in the familiar formula for the solution of a quadratic equation, x= −b + √ b2 − 4ac 2a (1.3.2) b 2 . (In §5.6 we

the addition becomes delicate and roundoff-prone whenever ac will learn how to avoid the problem in this particular case.)

30

Chapter 1.

Preliminaries

Roundoff error is a characteristic of computer hardware. There is another, different, kind of error that is a characteristic of the program or algorithm used, independent of the hardware on which the program is executed. Many numerical algorithms compute “discrete” approximations to some desired “continuous” quantity. For example, an integral is evaluated numerically by computing a function at a discrete set of points, rather than at “every” point. Or, a function may be evaluated by summing a ﬁnite number of leading terms in its inﬁnite series, rather than all inﬁnity terms. In cases like this, there is an adjustable parameter, e.g., the number of points or of terms, such that the “true” answer is obtained only when that parameter goes to inﬁnity. Any practical calculation is done with a ﬁnite, but sufﬁciently large, choice of that parameter. The discrepancy between the true answer and the answer obtained in a practical calculation is called the truncation error. Truncation error would persist even on a hypothetical, “perfect” computer that had an inﬁnitely accurate representation and no roundoff error. As a general rule there is not much that a programmer can do about roundoff error, other than to choose algorithms that do not magnify it unnecessarily (see discussion of “stability” below). Truncation error, on the other hand, is entirely under the programmer’s control. In fact, it is only a slight exaggeration to say that clever minimization of truncation error is practically the entire content of the ﬁeld of numerical analysis! Most of the time, truncation error and roundoff error do not strongly interact with one another. A calculation can be imagined as having, ﬁrst, the truncation error that it would have if run on an inﬁnite-precision computer, “plus” the roundoff error associated with the number of operations performed. Sometimes, however, an otherwise attractive method can be unstable. This means that any roundoff error that becomes “mixed into” the calculation at an early stage is successively magniﬁed until it comes to swamp the true answer. An unstable method would be useful on a hypothetical, perfect computer; but in this imperfect world it is necessary for us to require that algorithms be stable — or if unstable that we use them with great caution. Here is a simple, if somewhat artiﬁcial, example of an unstable algorithm: Suppose that it is desired to calculate all integer powers of the so-called “Golden Mean,” the number given by √ 5−1 ≈ 0.61803398 (1.3.3) φ≡ 2 It turns out (you can easily verify) that the powers φ n satisfy a simple recursion relation, φn+1 = φn−1 − φn (1.3.4)

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 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).

Thus, knowing the ﬁrst two values φ 0 = 1 and φ1 = 0.61803398, we can successively apply (1.3.4) performing only a single subtraction, rather than a slower multiplication by φ, at each stage. Unfortunately, the recurrence (1.3.4) also has another solution, namely the value √ − 1 ( 5 + 1). Since the recurrence is linear, and since this undesired solution has 2 magnitude greater than unity, any small admixture of it introduced by roundoff errors will grow exponentially. On a typical machine with 32-bit wordlength, (1.3.4) starts

1.3 Error, Accuracy, and Stability

31

to give completely wrong answers by about n = 16, at which point φ n is down to only 10−4 . The recurrence (1.3.4) is unstable, and cannot be used for the purpose stated. We will encounter the question of stability in many more sophisticated guises, later in this book.
CITED REFERENCES AND FURTHER READING: Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), Chapter 1. Kahaner, D., Moler, C., and Nash, S. 1989, Numerical Methods and Software (Englewood Cliffs, NJ: Prentice Hall), Chapter 2. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), §1.3. Wilkinson, J.H. 1964, Rounding Errors in Algebraic Processes (Englewood Cliffs, NJ: PrenticeHall).
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 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).

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 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).

Chapter 2. Solution of Linear Algebraic Equations

2.0 Introduction
A set of linear algebraic equations looks like this: a11 x1 + a12 x2 + a13 x3 + · · · + a1N xN = b1 a21 x1 + a22 x2 + a23 x3 + · · · + a2N xN = b2 a31 x1 + a32 x2 + a33 x3 + · · · + a3N xN = b3 ··· ··· aM1 x1 + aM2 x2 + aM3 x3 + · · · + aMN xN = bM Here the N unknowns x j , j = 1, 2, . . . , N are related by M equations. The coefﬁcients aij with i = 1, 2, . . . , M and j = 1, 2, . . . , N are known numbers, as are the right-hand side quantities b i , i = 1, 2, . . . , M . (2.0.1)

Nonsingular versus Singular Sets of Equations
If N = M then there are as many equations as unknowns, and there is a good chance of solving for a unique solution set of x j ’s. Analytically, there can fail to be a unique solution if one or more of the M equations is a linear combination of the others, a condition called row degeneracy, or if all equations contain certain variables only in exactly the same linear combination, called column degeneracy. (For square matrices, a row degeneracy implies a column degeneracy, and vice versa.) A set of equations that is degenerate is called singular. We will consider singular matrices in some detail in §2.6. Numerically, at least two additional things can go wrong: • While not exact linear combinations of each other, some of the equations may be so close to linearly dependent that roundoff errors in the machine render them linearly dependent at some stage in the solution process. In this case your numerical procedure will fail, and it can tell you that it has failed. 32

2.0 Introduction

33

• Accumulated roundoff errors in the solution process can swamp the true solution. This problem particularly emerges if N is too large. The numerical procedure does not fail algorithmically. However, it returns a set of x’s that are wrong, as can be discovered by direct substitution back into the original equations. The closer a set of equations is to being singular, the more likely this is to happen, since increasingly close cancellations will occur during the solution. In fact, the preceding item can be viewed as the special case where the loss of signiﬁcance is unfortunately total. Much of the sophistication of complicated “linear equation-solving packages” is devoted to the detection and/or correction of these two pathologies. As you work with large linear sets of equations, you will develop a feeling for when such sophistication is needed. It is difﬁcult to give any ﬁrm guidelines, since there is no such thing as a “typical” linear problem. But here is a rough idea: Linear sets with N as large as 20 or 50 can be routinely solved in single precision (32 bit ﬂoating representations) without resorting to sophisticated methods, if the equations are not close to singular. With double precision (60 or 64 bits), this number can readily be extended to N as large as several hundred, after which point the limiting factor is generally machine time, not accuracy. Even larger linear sets, N in the thousands or greater, can be solved when the coefﬁcients are sparse (that is, mostly zero), by methods that take advantage of the sparseness. We discuss this further in §2.7. At the other end of the spectrum, one seems just as often to encounter linear problems which, by their underlying nature, are close to singular. In this case, you might need to resort to sophisticated methods even for the case of N = 10 (though rarely for N = 5). Singular value decomposition (§2.6) is a technique that can sometimes turn singular problems into nonsingular ones, in which case additional sophistication becomes unnecessary.

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 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
Equation (2.0.1) can be written in matrix form as A·x=b (2.0.2)

Here the raised dot denotes matrix multiplication, A is the matrix of coefﬁcients, and b is the right-hand side written as a column vector, a11  a21 A= aM1  a12 a22 ··· aM2 ... ...  a1N a2N    b1 b  b= 2  ··· bM 

(2.0.3)

. . . aMN

By convention, the ﬁrst index on an element a ij denotes its row, the second index its column. For most purposes you don’t need to know how a matrix is stored in a computer’s physical memory; you simply reference matrix elements by their two-dimensional addresses, e.g., a 34 = a[3][4]. We have already seen, in §1.2, that this C notation can in fact hide a rather subtle and versatile physical storage scheme, “pointer to array of pointers to rows.” You might wish to review that section

34

Chapter 2.

Solution of Linear Algebraic Equations

at this point. Occasionally it is useful to be able to peer through the veil, for example to pass a whole row a[i][j], j=1, . . . , N by the reference a[i].

Tasks of Computational Linear Algebra
We will consider the following tasks as falling in the general purview of this chapter: • Solution of the matrix equation A·x = b for an unknown vector x, where A is a square matrix of coefﬁcients, raised dot denotes matrix multiplication, and b is a known right-hand side vector (§2.1–§2.10). • Solution of more than one matrix equation A · x j = bj , for a set of vectors xj , j = 1, 2, . . . , each corresponding to a different, known right-hand side vector bj . In this task the key simpliﬁcation is that the matrix A is held constant, while the right-hand sides, the b’s, are changed (§2.1–§2.10). • Calculation of the matrix A −1 which is the matrix inverse of a square matrix A, i.e., A · A−1 = A−1 · A = 1, where 1 is the identity matrix (all zeros except for ones on the diagonal). This task is equivalent, for an N × N matrix A, to the previous task with N different b j ’s (j = 1, 2, . . . , N ), namely the unit vectors (b j = all zero elements except for 1 in the jth component). The corresponding x’s are then the columns of the matrix inverse of A (§2.1 and §2.3). • Calculation of the determinant of a square matrix A (§2.3). If M < N , or if M = N but the equations are degenerate, then there are effectively fewer equations than unknowns. In this case there can be either no solution, or else more than one solution vector x. In the latter event, the solution space consists of a particular solution x p added to any linear combination of (typically) N − M vectors (which are said to be in the nullspace of the matrix A). The task of ﬁnding the solution space of A involves • Singular value decomposition of a matrix A. This subject is treated in §2.6. In the opposite case there are more equations than unknowns, M > N . When this occurs there is, in general, no solution vector x to equation (2.0.1), and the set of equations is said to be overdetermined. It happens frequently, however, that the best “compromise” solution is sought, the one that comes closest to satisfying all equations simultaneously. If closeness is deﬁned in the least-squares sense, i.e., that the sum of the squares of the differences between the left- and right-hand sides of equation (2.0.1) be minimized, then the overdetermined linear problem reduces to a (usually) solvable linear problem, called the • Linear least-squares problem. The reduced set of equations to be solved can be written as the N ×N set of equations (AT · A) · x = (AT · b) (2.0.4)
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 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).

where AT denotes the transpose of the matrix A. Equations (2.0.4) are called the normal equations of the linear least-squares problem. There is a close connection

2.0 Introduction

35

between singular value decomposition and the linear least-squares problem, and the latter is also discussed in §2.6. You should be warned that direct solution of the normal equations (2.0.4) is not generally the best way to ﬁnd least-squares solutions. Some other topics in this chapter include • Iterative improvement of a solution (§2.5) • Various special forms: symmetric positive-deﬁnite (§2.9), tridiagonal (§2.4), band diagonal (§2.4), Toeplitz (§2.8), Vandermonde (§2.8), sparse (§2.7) • Strassen’s “fast matrix inversion” (§2.11).
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 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).

Standard Subroutine Packages
We cannot hope, in this chapter or in this book, to tell you everything there is to know about the tasks that have been deﬁned above. In many cases you will have no alternative but to use sophisticated black-box program packages. Several good ones are available, though not always in C. LINPACK was developed at Argonne National Laboratories and deserves particular mention because it is published, documented, and available for free use. A successor to LINPACK, LAPACK, is now becoming available. Packages available commercially (though not necessarily in C) include those in the IMSL and NAG libraries. You should keep in mind that the sophisticated packages are designed with very large linear systems in mind. They therefore go to great effort to minimize not only the number of operations, but also the required storage. Routines for the various tasks are usually provided in several versions, corresponding to several possible simpliﬁcations in the form of the input coefﬁcient matrix: symmetric, triangular, banded, positive deﬁnite, etc. If you have a large matrix in one of these forms, you should certainly take advantage of the increased efﬁciency provided by these different routines, and not just use the form provided for general matrices. There is also a great watershed dividing routines that are direct (i.e., execute in a predictable number of operations) from routines that are iterative (i.e., attempt to converge to the desired answer in however many steps are necessary). Iterative methods become preferable when the battle against loss of signiﬁcance is in danger of being lost, either due to large N or because the problem is close to singular. We will treat iterative methods only incompletely in this book, in §2.7 and in Chapters 18 and 19. These methods are important, but mostly beyond our scope. We will, however, discuss in detail a technique which is on the borderline between direct and iterative methods, namely the iterative improvement of a solution that has been obtained by direct methods (§2.5).
CITED REFERENCES AND FURTHER READING: Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press). Gill, P.E., Murray, W., and Wright, M.H. 1991, Numerical Linear Algebra and Optimization, vol. 1 (Redwood City, CA: Addison-Wesley). Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), Chapter 4. Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.).

36

Chapter 2.

Solution of Linear Algebraic Equations

Coleman, T.F., and Van Loan, C. 1988, Handbook for Matrix Computations (Philadelphia: S.I.A.M.). Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Englewood Cliffs, NJ: Prentice-Hall). Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations (New York: Wiley). Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), Chapter 2. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York: McGraw-Hill), Chapter 9.

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 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).

2.1 Gauss-Jordan Elimination
For inverting a matrix, Gauss-Jordan elimination is about as efﬁcient as any other method. For solving sets of linear equations, Gauss-Jordan elimination produces both the solution of the equations for one or more right-hand side vectors b, and also the matrix inverse A −1 . However, its principal weaknesses are (i) that it requires all the right-hand sides to be stored and manipulated at the same time, and (ii) that when the inverse matrix is not desired, Gauss-Jordan is three times slower than the best alternative technique for solving a single linear set (§2.3). The method’s principal strength is that it is as stable as any other direct method, perhaps even a bit more stable when full pivoting is used (see below). If you come along later with an additional right-hand side vector, you can multiply it by the inverse matrix, of course. This does give an answer, but one that is quite susceptible to roundoff error, not nearly as good as if the new vector had been included with the set of right-hand side vectors in the ﬁrst instance. For these reasons, Gauss-Jordan elimination should usually not be your method of ﬁrst choice, either for solving linear equations or for matrix inversion. The decomposition methods in §2.3 are better. Why do we give you Gauss-Jordan at all? Because it is straightforward, understandable, solid as a rock, and an exceptionally good “psychological” backup for those times that something is going wrong and you think it might be your linear-equation solver. Some people believe that the backup is more than psychological, that GaussJordan elimination is an “independent” numerical method. This turns out to be mostly myth. Except for the relatively minor differences in pivoting, described below, the actual sequence of operations performed in Gauss-Jordan elimination is very closely related to that performed by the routines in the next two sections. For clarity, and to avoid writing endless ellipses (· · ·) we will write out equations only for the case of four equations and four unknowns, and with three different righthand side vectors that are known in advance. You can write bigger matrices and extend the equations to the case of N × N matrices, with M sets of right-hand side vectors, in completely analogous fashion. The routine implemented below is, of course, general.

2.1 Gauss-Jordan Elimination

37

Elimination on Column-Augmented Matrices

a11  a21 a31 a41

Consider the linear matrix equation      a12 a22 a32 a42 a13 a23 a33 a43 a14 x11 a24   x21  · a34 x31 a44 x41 x12  x22  x32 x42



x13  x23  x33 x43





y11  y21 y31 y41 0 1 0 0 0 0 1 0

y12 y22 y32 y42

y13 y23 y33 y43

y14 y24  y34 y44


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 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).

 =

b11  b21  b31 b41





b12  b22  b32 b42





b13  b23  b33 b43





1 0 0 0

0 0  0 1

 (2.1.1)

Here the raised dot (·) signiﬁes matrix multiplication, while the operator just signiﬁes column augmentation, that is, removing the abutting parentheses and making a wider matrix out of the operands of the operator. It should not take you long to write out equation (2.1.1) and to see that it simply states that xij is the ith component (i = 1, 2, 3, 4) of the vector solution of the jth right-hand side (j = 1, 2, 3), the one whose coefﬁcients are b ij , i = 1, 2, 3, 4; and that the matrix of unknown coefﬁcients y ij is the inverse matrix of a ij . In other words, the matrix solution of [A] · [x1 x2 x3 Y] = [b1 b2 b3 1] (2.1.2)

where A and Y are square matrices, the b i ’s and xi ’s are column vectors, and 1 is the identity matrix, simultaneously solves the linear sets A · x1 = b1 and A·Y=1 (2.1.4) Now it is also elementary to verify the following facts about (2.1.1): • Interchanging any two rows of A and the corresponding rows of the b’s and of 1, does not change (or scramble in any way) the solution x’s and Y. Rather, it just corresponds to writing the same set of linear equations in a different order. • Likewise, the solution set is unchanged and in no way scrambled if we replace any row in A by a linear combination of itself and any other row, as long as we do the same linear combination of the rows of the b’s and 1 (which then is no longer the identity matrix, of course). • Interchanging any two columns of A gives the same solution set only if we simultaneously interchange corresponding rows of the x’s and of Y. In other words, this interchange scrambles the order of the rows in the solution. If we do this, we will need to unscramble the solution by restoring the rows to their original order. Gauss-Jordan elimination uses one or more of the above operations to reduce the matrix A to the identity matrix. When this is accomplished, the right-hand side becomes the solution set, as one sees instantly from (2.1.2). A · x2 = b2 A · x3 = b3 (2.1.3)

38

Chapter 2.

Solution of Linear Algebraic Equations

Pivoting
In “Gauss-Jordan elimination with no pivoting,” only the second operation in the above list is used. The ﬁrst row is divided by the element a 11 (this being a trivial linear combination of the ﬁrst row with any other row — zero coefﬁcient for the other row). Then the right amount of the ﬁrst row is subtracted from each other row to make all the remaining a i1 ’s zero. The ﬁrst column of A now agrees with the identity matrix. We move to the second column and divide the second row by a22 , then subtract the right amount of the second row from rows 1, 3, and 4, so as to make their entries in the second column zero. The second column is now reduced to the identity form. And so on for the third and fourth columns. As we do these operations to A, we of course also do the corresponding operations to the b’s and to 1 (which by now no longer resembles the identity matrix in any way!). Obviously we will run into trouble if we ever encounter a zero element on the (then current) diagonal when we are going to divide by the diagonal element. (The element that we divide by, incidentally, is called the pivot element or pivot.) Not so obvious, but true, is the fact that Gauss-Jordan elimination with no pivoting (no use of the ﬁrst or third procedures in the above list) is numerically unstable in the presence of any roundoff error, even when a zero pivot is not encountered. You must never do Gauss-Jordan elimination (or Gaussian elimination, see below) without pivoting! So what is this magic pivoting? Nothing more than interchanging rows (partial pivoting) or rows and columns (full pivoting), so as to put a particularly desirable element in the diagonal position from which the pivot is about to be selected. Since we don’t want to mess up the part of the identity matrix that we have already built up, we can choose among elements that are both (i) on rows below (or on) the one that is about to be normalized, and also (ii) on columns to the right (or on) the column we are about to eliminate. Partial pivoting is easier than full pivoting, because we don’t have to keep track of the permutation of the solution vector. Partial pivoting makes available as pivots only the elements already in the correct column. It turns out that partial pivoting is “almost” as good as full pivoting, in a sense that can be made mathematically precise, but which need not concern us here (for discussion and references, see [1]). To show you both variants, we do full pivoting in the routine in this section, partial pivoting in §2.3. We have to state how to recognize a particularly desirable pivot when we see one. The answer to this is not completely known theoretically. It is known, both theoretically and in practice, that simply picking the largest (in magnitude) available element as the pivot is a very good choice. A curiosity of this procedure, however, is that the choice of pivot will depend on the original scaling of the equations. If we take the third linear equation in our original set and multiply it by a factor of a million, it is almost guaranteed that it will contribute the ﬁrst pivot; yet the underlying solution of the equations is not changed by this multiplication! One therefore sometimes sees routines which choose as pivot that element which would have been largest if the original equations had all been scaled to have their largest coefﬁcient normalized to unity. This is called implicit pivoting. There is some extra bookkeeping to keep track of the scale factors by which the rows would have been multiplied. (The routines in §2.3 include implicit pivoting, but the routine in this section does not.) Finally, let us consider the storage requirements of the method. With a little reﬂection you will see that at every stage of the algorithm, either an element of A is

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 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).

2.1 Gauss-Jordan Elimination

39

predictably a one or zero (if it is already in a part of the matrix that has been reduced to identity form) or else the exactly corresponding element of the matrix that started as 1 is predictably a one or zero (if its mate in A has not been reduced to the identity form). Therefore the matrix 1 does not have to exist as separate storage: The matrix inverse of A is gradually built up in A as the original A is destroyed. Likewise, the solution vectors x can gradually replace the right-hand side vectors b and share the same storage, since after each column in A is reduced, the corresponding row entry in the b’s is never again used. Here is the routine for Gauss-Jordan elimination with full pivoting:
#include #include "nrutil.h" #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} void gaussj(float **a, int n, float **b, int m) Linear equation solution by Gauss-Jordan elimination, equation (2.1.1) above. a[1..n][1..n] is the input matrix. b[1..n][1..m] is input containing the m right-hand side vectors. On output, a is replaced by its matrix inverse, and b is replaced by the corresponding set of solution vectors. { int *indxc,*indxr,*ipiv; int i,icol,irow,j,k,l,ll; float big,dum,pivinv,temp; indxc=ivector(1,n); The integer arrays ipiv, indxr, and indxc are indxr=ivector(1,n); used for bookkeeping on the pivoting. ipiv=ivector(1,n); for (j=1;j

Numerical Recipes in C

Similar Documents

Linear Progr

Cholesky Decomposition

Managerial Economics

Thesis Referrence

Investigatory Project

Python Programming

Programming Languages

Arduino

Asian Options

Restaurants

Eigen Value Eigen Vectors

Questionnaire on Mobile Operator

Dates Palm Research Paper

Swag

Unions

Popular Essays