f0-0

Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Numerical Recipes in Fortran 77 The Art of ScientificComputing Second Edition Volume 1 of Fortran NumericalRecipes WilliamH.Press Harvard-SmithsonianCenterforAstrophysics SaulA.Teukolsky DepartmentofPhysics,CornellUniversity WilliamT. Vetterling PolaroidCorporation Brian P. Flannery EXXONResearchandEngineeringCompany Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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, USA10 Stamford Road, Oakleigh, Melbourne 3166, Australia Copyright c/circlecopyrtCambridge University Press 1986, 1992 except for §13.10, which is placed into the public domain, and except for all other computer programs and procedures, which are Copyright c/circlecopyrtNumerical Recipes Software 1986, 1992, 1997 All Rights Reserved. Some sections of this book were originally published, in different form, in Computers in Physics magazine, Copyright c/circlecopyrtAmerican Institute of Physics, 1988–1992. First Edition originally published 1986; Second Edition originally published 1992 as Numerical Recipes in FORTRAN: The Art of Scientific Computing Reprinted with corrections, 1993, 1994, 1995. Reprinted with corrections, 1996, 1997, 2001, as Numerical Recipes in Fortran 77: The Art of Scientific Computing (Vol. 1 of Fortran Numerical Recipes) This reprinting is corrected to software version 2.10 Printed in the United States of America Typeset in T EX 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 thesoftware 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. xxi. See the section “License Information” (pp.xx–xxiii) forinformation onobtaining moregeneral licenses atlow cost. Machine-readable mediacontaining thesoftware inthisbook,withincluded licenses for use on a single screen, are available from Cambridge University Press. See theorder form at the back of the book, email to “[email protected]” (North America) or“[email protected]” (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 computerexceptonethatisspecificallylicensed, isstrictlyprohibited. Technicalquestions,corrections, and requests for information should be addressed to Numerical RecipesSoftware, P.O. Box 243, Cambridge, MA 02238 (USA), email “[email protected]”, or fax 781 863-1739. Library of Congress Cataloging in Publication Data Numerical recipes in Fortran 77 : the art of scientific computing / William H. Press ...[et al.]. – 2nd ed. Includes bibliographical references (p. ) and index.ISBN 0-521-43064-X 1. Numerical analysis–Computer programs. 2. Science–Mathematics–Computer programs. 3. FORTRAN (Computer program language) I. Press, William H. QA297.N866 1992519.4 /prime0285/prime53–dc20 92-8876 A catalog record for this book is available from the British Library. ISBN 0 521 43064 X Volume 1 (this book) ISBN 0 521 57439 0 Volume 2 ISBN 0 521 43721 0 Example book in FORTRAN ISBN 0 521 57440 4 FORTRAN diskette (IBM 3.5/prime/prime) ISBN 0 521 57608 3 CDROM (IBM PC/Macintosh) ISBN 0 521 57607 5 CDROM (UNIX) Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Contents Plan of the Two-VolumeEdition xiii Preface tothe Second Edition xv Preface tothe FirstEdition xviiiLicense Information xx ComputerPrograms by Chapter and Section xxiv 1 Preliminaries 1 1.0Introduction 1 1.1ProgramOrganizationandControlStructures 5 1.2Error,Accuracy,andStability 18 2 Solutionof Linear AlgebraicEquations 22 2.0Introduction 222.1Gauss-JordanElimination 27 2.2GaussianEliminationwithBacksubstitution 33 2.3LUDecompositionandIts Applications 34 2.4TridiagonalandBandDiagonalSystemsofEquations 422.5IterativeImprovementofa SolutiontoLinearEquations 47 2.6SingularValueDecomposition 51 2.7SparseLinearSystems 632.8VandermondeMatrices andToeplitzMatrices 82 2.9CholeskyDecomposition 89 2.10QR Decomposition 91 2.11Is MatrixInversionan N 3Process? 95 3 Interpolationand Extrapolation 99 3.0Introduction 99 3.1PolynomialInterpolationandExtrapolation 102 3.2RationalFunctionInterpolationandExtrapolation 104 3.3CubicSplineInterpolation 107 3.4HowtoSearchanOrderedTable 1103.5Coefficientsofthe InterpolatingPolynomial 113 3.6Interpolationin TwoorMoreDimensions 116 v vi ContentsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).4 Integration ofFunctions 123 4.0Introduction 123 4.1Classical FormulasforEquallySpacedAbscissas 1244.2ElementaryAlgorithms 130 4.3RombergIntegration 134 4.4ImproperIntegrals 135 4.5GaussianQuadraturesandOrthogonalPolynomials 140 4.6MultidimensionalIntegrals 155 5 Evaluation ofFunctions 159 5.0Introduction 159 5.1SeriesandTheirConvergence 159 5.2EvaluationofContinuedFractions 163 5.3PolynomialsandRationalFunctions 167 5.4ComplexArithmetic 1715.5RecurrenceRelations andClenshaw’s RecurrenceFormula 172 5.6QuadraticandCubicEquations 178 5.7NumericalDerivatives 1805.8ChebyshevApproximation 184 5.9DerivativesorIntegralsofa Chebyshev-approximatedFunction 189 5.10PolynomialApproximationfromChebyshevCoefficients 191 5.11EconomizationofPowerSeries 192 5.12Pad´eApproximants 194 5.13RationalChebyshevApproximation 197 5.14EvaluationofFunctionsbyPath Integration 201 6 Special Functions 205 6.0Introduction 205 6.1GammaFunction,BetaFunction,Factorials,BinomialCoefficients 206 6.2IncompleteGammaFunction,ErrorFunction,Chi-Square ProbabilityFunction,CumulativePoissonFunction 209 6.3ExponentialIntegrals 215 6.4IncompleteBeta Function,Student’sDistribution,F-Distribution, CumulativeBinomialDistribution 219 6.5Bessel FunctionsofIntegerOrder 223 6.6ModifiedBessel FunctionsofIntegerOrder 2296.7Bessel FunctionsofFractionalOrder,AiryFunctions,Spherical Bessel Functions 234 6.8SphericalHarmonics 246 6.9FresnelIntegrals,CosineandSineIntegrals 248 6.10Dawson’sIntegral 2526.11EllipticIntegralsandJacobianEllipticFunctions 254 6.12HypergeometricFunctions 263 7 RandomNumbers 266 7.0Introduction 2667.1UniformDeviates 267 Contents viiSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).7.2TransformationMethod: ExponentialandNormalDeviates 277 7.3RejectionMethod: Gamma,Poisson,BinomialDeviates 2817.4GenerationofRandomBits 287 7.5RandomSequencesBased onData Encryption 290 7.6SimpleMonteCarloIntegration 2957.7Quasi-(thatis,Sub-)RandomSequences 299 7.8AdaptiveandRecursiveMonteCarloMethods 306 8 Sorting 320 8.0Introduction 3208.1StraightInsertionandShell’s Method 321 8.2Quicksort 3238.3Heapsort 327 8.4IndexingandRanking 329 8.5Selectingthe MthLargest 333 8.6DeterminationofEquivalenceClasses 337 9 Root Findingand NonlinearSets of Equations 340 9.0Introduction 3409.1BracketingandBisection 343 9.2SecantMethod,False PositionMethod,andRidders’Method 3479.3VanWijngaarden–Dekker–BrentMethod 352 9.4Newton-RaphsonMethodUsingDerivative 355 9.5RootsofPolynomials 3629.6Newton-RaphsonMethodforNonlinearSystemsofEquations 372 9.7GloballyConvergentMethodsforNonlinearSystemsofEquations 376 10 MinimizationorMaximization ofFunctions 387 10.0Introduction 38710.1GoldenSectionSearchinOneDimension 390 10.2ParabolicInterpolationandBrent’s MethodinOne Dimension 39510.3One-DimensionalSearchwithFirst Derivatives 399 10.4DownhillSimplexMethodinMultidimensions 402 10.5DirectionSet (Powell’s)MethodsinMultidimensions 406 10.6ConjugateGradientMethodsin Multidimensions 413 10.7VariableMetric MethodsinMultidimensions 41810.8LinearProgrammingandtheSimplexMethod 423 10.9SimulatedAnnealingMethods 436 11 Eigensystems 449 11.0Introduction 44911.1JacobiTransformationsofaSymmetricMatrix 456 11.2Reductionofa SymmetricMatrixtoTridiagonalForm: GivensandHouseholderReductions 462 11.3EigenvaluesandEigenvectorsofa TridiagonalMatrix 469 11.4HermitianMatrices 475 11.5Reductionofa GeneralMatrixtoHessenbergForm 476 viii ContentsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).11.6TheQR AlgorithmforReal HessenbergMatrices 480 11.7ImprovingEigenvaluesand/orFindingEigenvectorsby InverseIteration 487 12 Fast FourierTransform 490 12.0Introduction 49012.1FourierTransformofDiscretelySampledData 494 12.2Fast FourierTransform(FFT) 498 12.3FFT ofReal Functions,SineandCosineTransforms 50412.4FFT inTwo orMoreDimensions 515 12.5FourierTransformsofRealData inTwo andThreeDimensions 519 12.6ExternalStorageorMemory-LocalFFTs 525 13 Fourierand Spectral Applications 530 13.0Introduction 53013.1ConvolutionandDeconvolutionUsingtheFFT 53113.2CorrelationandAutocorrelationUsingtheFFT 538 13.3Optimal(Wiener)FilteringwiththeFFT 539 13.4PowerSpectrumEstimationUsing theFFT 542 13.5DigitalFilteringin theTimeDomain 551 13.6LinearPredictionandLinearPredictiveCoding 55713.7PowerSpectrumEstimationbytheMaximumEntropy (All Poles)Method 565 13.8SpectralAnalysisofUnevenlySampledData 56913.9ComputingFourierIntegralsUsing theFFT 577 13.10WaveletTransforms 584 13.11NumericalUse oftheSamplingTheorem 600 14 Statistical Description of Data 603 14.0Introduction 60314.1Momentsofa Distribution: Mean,Variance,Skewness, andSoForth 604 14.2DoTwo DistributionsHave theSameMeansorVariances? 609 14.3AreTwo DistributionsDifferent? 614 14.4ContingencyTableAnalysisofTwo Distributions 622 14.5LinearCorrelation 63014.6NonparametricorRankCorrelation 633 14.7DoTwo-DimensionalDistributionsDiffer? 640 14.8Savitzky-GolaySmoothingFilters 644 15 Modelingof Data 650 15.0Introduction 65015.1Least Squaresas a MaximumLikelihoodEstimator 65115.2FittingData toa StraightLine 655 15.3Straight-LineDatawithErrorsinBothCoordinates 660 15.4GeneralLinearLeast Squares 665 15.5NonlinearModels 675 Contents ixSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).15.6ConfidenceLimits onEstimatedModelParameters 684 15.7RobustEstimation 694 16 Integration ofOrdinaryDifferential Equations 701 16.0Introduction 701 16.1Runge-KuttaMethod 704 16.2AdaptiveStepsizeControlforRunge-Kutta 70816.3ModifiedMidpointMethod 716 16.4RichardsonExtrapolationandthe Bulirsch-StoerMethod 718 16.5Second-OrderConservativeEquations 72616.6StiffSets ofEquations 727 16.7Multistep,Multivalue,andPredictor-CorrectorMethods 740 17 TwoPoint BoundaryValue Problems 745 17.0Introduction 74517.1TheShootingMethod 749 17.2Shootingtoa FittingPoint 75117.3RelaxationMethods 753 17.4A WorkedExample: SpheroidalHarmonics 764 17.5AutomatedAllocationofMeshPoints 77417.6HandlingInternalBoundaryConditionsorSingularPoints 775 18 Integral Equations and Inverse Theory 779 18.0Introduction 77918.1FredholmEquationsofthe SecondKind 782 18.2VolterraEquations 786 18.3IntegralEquationswithSingularKernels 78818.4InverseProblemsandtheUse ofAPrioriInformation 795 18.5LinearRegularizationMethods 799 18.6Backus-GilbertMethod 80618.7MaximumEntropyImageRestoration 809 19 Partial Differential Equations 818 19.0Introduction 81819.1Flux-ConservativeInitial ValueProblems 825 19.2DiffusiveInitialValueProblems 838 19.3InitialValueProblemsin Multidimensions 84419.4FourierandCyclicReductionMethodsforBoundary ValueProblems 848 19.5RelaxationMethodsforBoundaryValueProblems 85419.6MultigridMethodsforBoundaryValueProblems 862 20 Less-Numerical Algorithms 881 20.0Introduction 88120.1DiagnosingMachineParameters 881 20.2GrayCodes 886 x ContentsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).20.3Cyclic RedundancyandOtherChecksums 888 20.4HuffmanCodingandCompressionofData 89620.5ArithmeticCoding 902 20.6ArithmeticatArbitraryPrecision 906 References forVolume 1 916 Indexof Programs and Dependencies (Vol. 1) 921 General Indexto Volumes 1 and 2 Contents of Volume 2: NumericalRecipes inFortran 90 Preface toVolume 2 viii Forewordby Michael Metcalf x License Information xvii 21 Introductionto Fortran90 Language Features 935 22 Introductionto Parallel Programming 96223 NumericalRecipes Utilities forFortran 90 987 Fortran90 Code Chapters 1009 B1 Preliminaries 1010 B2 Solutionof Linear AlgebraicEquations 1014B3 Interpolationand Extrapolation 1043 B4 Integration ofFunctions 1052 B5 Evaluation ofFunctions 1070B6 Special Functions 1083B7 RandomNumbers 1141 B8 Sorting 1167 B9 Root Findingand NonlinearSets of Equations 1182B10 MinimizationorMaximization ofFunctions 1201B11 Eigensystems 1225 B12 Fast FourierTransform 1235 Contents xiSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).B13 Fourierand Spectral Applications 1253 B14 Statistical Description of Data 1269B15 Modelingof Data 1285B16 Integration ofOrdinaryDifferential Equations 1297 B17 TwoPoint BoundaryValue Problems 1314 B18 Integral Equations and Inverse Theory 1325B19 Partial Differential Equations 1332B20 Less-Numerical Algorithms 1343 References forVolume 2 1359Appendices C1 Listingof Utility Modules(nrtype andnrutil) 1361 C2 Listingof Explicit Interfaces 1384C3 Indexof Programs and Dependencies (Vol. 2) 1434 General Indexto Volumes 1 and 2 1447 Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica). xii Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Planofthe Two-Volume Edition Fortran,longtheepitomeofstability,isonceagainalanguageinflux. Fortran90 is not just the long-awaited updatingof traditional Fortran 77 to moderncomputing practices, but also demonstrates Fortran’s decisive bid to be the language of choice for parallel programming on multiprocessor computers. At the same time, Fortran 90 is completely backwards-compatible with all Fortran 77 code. So, users with legacy code, or who choose to use only older language constructs, will still get the benefit of updated and actively maintained compilers. As we, the authors of Numerical Recipes , watched the gestation and birth of Fortran 90 by its governing standards committee (an interesting process described by a leading Committee member, Michael Metcalf, in the Foreword to our Volume 2), it became clear to us that the right moment for movingNumerical Recipes fromFortran 77 to Fortran 90 was sooner, rather than later. On the other hand, it was equally clear that Fortran-77-style programming — no matter whether with Fortran 77 or Fortran 90 compilers — is, and will continuefor a long time to be, the “mother tongue” of a large populationof active scientists, engineers, and other users of numerical computation. This is not a user base that we would willingly or knowingly abandon. The solution was immediately clear: a two-volume edition of the Fortran Numerical Recipes consisting of Volume 1 (this one, a corrected reprinting of the previous one-volume edition), now retitled Numerical Recipes in Fortran 77 , and a completelynewVolume2,titled NumericalRecipesinFortran90: TheArtof Parallel Scientific Computing . Volume 2 begins with three chapters (21, 22, and 23) that extend the narrative of the first volume to the new subjects of Fortran 90 language features, parallel programming methodology, and the implementation of certain useful utility functions in Fortran 90. Then, in exact correspondence with Volume 1’s Chapters 1–20, are new chapters B1–B20, devoted principally to the listing and explanation of new Fortran 90 routines. With a few exceptions, each Fortran 77routinein Volume 1 has a correspondingnew Fortran 90 versionin Volume 2. (The exceptions are a few new capabilities, notably in random number generation and in multigrid PDE solvers, that are unique to Volume 2’s Fortran 90.) Otherwise, thereis no duplication between the volumes. The detailed explanation of the algorithms in this Volume 1 is intended to apply to, and be essential for, both volumes. Inotherwords: YoucanusethisVolume1withouthavingVolume2,butyou can’tuseVolume2withoutVolume1. Wethinkthatthereismuchtobegainedby havingandusing bothvolumes: Fortran90’sparallellanguageconstructionsarenot only useful for present and future multiprocessor machines; they also allow for the elegant and concise formulation of many algorithms on ordinary single-processor computers. We think that essentially allFortran programmers will want gradually to migrate into Fortran 90 and into a mode of “thinking parallel.” We have written Volume 2 specifically to help with this important transition. Volume 2’s discussion of parallel programming is focused on those issues of directrelevancetotheFortran90programmer. Somemoregeneralaspectsofparallel programming,suchascommunicationcosts,synchronizationofmultipleprocessers,etc., are touched on only briefly. We provide references to the extensive literature on these more specialized topics. xiii xiv PlanoftheTwo-VolumeEditionSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).A special note to C programmers: Right now, there is no effort at producing a parallel version of C that is comparable to Fortran 90 in maturity, acceptance,and stability. We think, therefore, that C programmers will be well served by using Volume 2, either in conjuction with this Volume 1 or else in conjunctionwith the sister volume Numerical Recipes in C: The Art of Scientific Computing , for an educational excursion into Fortran 90, its parallel programming constructions, and the numerical algorithms that capitalize on them. C and C++ programming have not been far from our minds as we have written this two-volume version. We think you will find that time spent in absorbing the principal lessons of Volume 2’s Chapters 21–23 will be amply repaid in the future, as C and C++ eventuallydevelop standard parallel extensions. Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Prefaceto 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 first edition puts us now in a difficult,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 dangerthat,ifwearenotcareful,wemightproduceasecondeditionthatisweighty, balanced, scholarly, and boring. Itis amixedblessingthatweknowmorenowthanwedidsixyearsago. Then, weweremakingeducatedguesses,basedonexistingliteratureandourownresearch, aboutwhichnumericaltechniqueswerethemostimportantandrobust. Now,wehave thebenefitofdirectfeedbackfromalargereadercommunity. Letterstoouralter-egoenterprise,NumericalRecipesSoftware,areinthethousandsperyear. (Please, don’t telephone us.) Our post office box has become a magnet for letters pointing out that we have omitted some particular technique, well known to be important in aparticular field of science or engineering. We value such letters, and digest them carefully,especially when they point us to specific referencesin the literature. The inevitable result of this input is that this Second Edition of Numerical Recipesis substantially larger than its predecessor, in fact about 50% larger both in words and number of included programs (the latter now numberingwell 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 increasedin difficulty,the book’sprincipal 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, find 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´e approximants, and rational Chebyshev approximation Bessel functions, and modified 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 high- dimensional spaces globally convergent methods for sets of nonlinear equations xv xvi PrefacetotheSecondEditionSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).simulated annealing minimization for continuous control spaces fast Fouriertransform(FFT) forreal data in two andthree 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 filters fitting straight line data with errors in both coordinates a two-dimensional Kolmogorov-Smirnoff test the statistical bootstrap method embeddedRunge-Kutta-Fehlbergmethods for differential equations high-order methods for stiff differential equations a new chapter on “less-numerical” algorithms, including Huffman and arithmeticcoding,arbitraryprecisionarithmetic,andseveralothertopics. Consult the Preface to the First Edition, following, or the Table of Contents, for alist of the more “basic” subjects treated. Acknowledgments It is not possible for us to list by name here all the readers who have made usefulsuggestions; we are gratefulforthese. In the text, we attempt to givespecific 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,DouglasEardley,RameshNarayan,DavidSpergel,AlanOppenheim,SallieBaliunas, Scott Tremaine, Glennys Farrar, Steven Block, John Peacock, Thomas Loredo, Matthew Choptuik, Gregory Cook, L. Samuel Finn, P. Deuflhard, Harold Lewis, Peter Weinberger, David Syer, Richard Ferch, Steven Ebstein, 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 influenced many of the routines in this book, and wrote or rewrote many more routines in its C-languagetwin and the companionExample books. Our project has benefited enormously from Seth’s talent for detecting, and following the trail of, even very slight anomalies (often compiler bugs, but occasionally our errors), andfrom his good programming sense. We prepared this book for publication on DEC and Sun workstations run- ning the UNIX operating system, and on a 486/33 PC compatible runningMS-DOS5.0/Windows3.0. (See §1.0 for a list of additional computers used in program tests.) We enthusiastically recommend the principal software used: GNU Emacs, T EX, Perl, Adobe Illustrator, and PostScript. Also used were a variety ofFORTRAN compilers — too numerous (and sometimes too buggy) for individual PrefacetotheSecondEdition xviiSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).acknowledgment. It is a sobering fact that our standard test suite (exercisingall the routinesinthisbook)hasuncoveredcompilerbugsinalargemajorityofthecompil-ers tried. When possible, we work with developers to see that such bugs get fixed; weencourageinterestedcompilerdeveloperstocontactusaboutsucharrangements. WHPandSATacknowledgethecontinuedsupportoftheU.S.NationalScience 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.TeukolskyWilliamT.Vetterling BrianP. Flannery Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Preface tothe FirstEdition Wecallthisbook NumericalRecipes forseveralreasons. Inonesense,thisbook is indeed a “cookbook” on numerical computation. However there is an importantdistinction between a cookbook and a restaurant menu. The latter presents choices among complete dishes in each of which the individual flavors 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 imple-mentations of these ideas in the form of working computer routines. Our task has been to find the right balance among these ingredients for each topic. You will find that for some topics we have tilted quite far to the analytic side; this where wehave 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. Weadmit,therefore,tosomeunevennessinthe“level”ofthisbook. Abouthalf of it is suitable foran advancedundergraduatecourseon numericalcomputationfor science or engineering majors. The other half ranges from the level of a graduate courseto that ofa professionalreference. Most cookbookshave,afterall, recipes atvaryinglevelsofcomplexity. Anattractivefeatureofthisapproach,wethink,isthat thereadercanusethebookatincreasinglevelsofsophisticationashis/herexperience grows. Eveninexperiencedreadersshouldbeabletouseourmostadvancedroutines 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 efficient, clever, and — important — clear. The alternative viewpoint, that efficient computational methods mustnecessarily be so arcane and complex as to be useful only in “black box” form, we firmly 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 specific 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 assumethat 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, wedohave one introductory chapter on methods for partial differential equations (Chapter19). Second,weobviouslycannotinclude everything else. Alltheso-called “standard” topics of a numerical analysis course have been included in this book: xviii PrefacetotheFirst Edition xixSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).linear equations (Chapter 2), interpolationand extrapolation(Chaper 3), integration (Chaper 4), nonlinear root-finding (Chapter 9), eigensystems (Chapter 11), andordinary 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. Someothersubjectsthatwecoverindetailarenotusuallyfoundinthestandard numericalanalysis texts. These includethe evaluationoffunctionsandofparticular special functions of higher mathematics (Chapters 5 and 6); random numbers and Monte Carlo methods (Chapter 7); sorting (Chapter 8); optimization, including multidimensionalmethods (Chapter10); Fouriertransformmethods, includingFFTmethods 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 boundaryvalue problems, bothshooting and relaxationmethods (Chapter 17). TheprogramsinthisbookareincludedinANSI-standard FORTRAN-77. Versions of the book in C,Pascal, and BASICare available separately. We have more to say about the FORTRAN language, and the computational environment assumed by our routines, in §1.1 (Introduction). Acknowledgments Manycolleagueshavebeengenerousingivingusthebenefitoftheirnumerical and computational experience, in providing us with programs, in commenting on the manuscript,or in generalencouragement. We particularlywish to thank George Rybicki, Douglas Eardley, Philip Marcus, Stuart Shapiro, Paul Horowitz, BruceMusicus, Irwin Shapiro, Stephen Wolfram, Henry Abarbanel,Larry Smarr, Richard Muller, John Bahcall, and A.G.W. Cameron. Wealsowishtoacknowledgetwoindividualswhomwehavenevermet: Forman Acton,whose1970textbook Numerical MethodsthatWork (NewYork: Harperand Row) has surely left its stylistic mark on us; and Donald Knuth, both for his seriesof books on The Art of Computer Programming (Reading, MA: Addison-Wesley), and for TEX, the computertypesetting 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 BrianP. FlannerySaul A.Teukolsky WilliamT.Vetterling