Table Of ContentNumerical Recipes in C
The Art of Scientific Computing
Second Edition
WilliamH.Press
Harvard-SmithsonianCenterforAstrophysics
SaulA.Teukolsky
DepartmentofPhysics,CornellUniversity
WilliamT.Vetterling
PolaroidCorporation
BrianP.Flannery
EXXONResearchandEngineeringCompany
CAMBRIDGE UNIVERSITY PRESS
Cambridge NewYork PortChester Melbourne Sydney
PublishedbythePressSyndicateoftheUniversityofCambridge
ThePitt Building,TrumpingtonStreet, CambridgeCB21RP
40West 20thStreet, NewYork, NY10011-4211,USA
10StamfordRoad,Oakleigh,Melbourne3166,Australia
Copyright(cid:1)c CambridgeUniversityPress 1988,1992
exceptfor§13.10andAppendixB,whichareplacedintothepublicdomain,
andexceptforallothercomputerprogramsandprocedures,whichare
Copyright(cid:1)c NumericalRecipesSoftware1987,1988,1992,1997
All Rights Reserved.
Somesectionsofthisbookwereoriginallypublished,indifferentform,inComputers
inPhysicsmagazine,Copyright(cid:1)c AmericanInstituteofPhysics,1988–1992.
FirstEditionoriginallypublished1988;SecondEditionoriginallypublished1992.
Reprintedwith corrections,1993, 1994,1995,1997.
This reprintingis correctedtosoftwareversion2.08
Printed in the United States of America
Typeset in TEX
Withoutanadditionallicensetousethecontainedsoftware,thisbookisintendedas
atextandreferencebook,forreadingpurposesonly. Afreelicenseforlimiteduseofthe
softwarebytheindividualownerofacopyofthisbookwhopersonallytypesoneormore
routinesintoasinglecomputerisgrantedundertermsdescribedonp.xvii.Seethesection
“LicenseInformation”(pp.xvi–xviii)forinformationonobtainingmoregenerallicenses
atlowcost.
Machine-readablemediacontainingthesoftwareinthisbook,withincludedlicenses
for use on a single screen, are available from Cambridge University Press. See the
order 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
MidlandAvenue,PortChester,NY10573(USA),forfurtherinformation.
The software may also be downloaded, with immediate purchase of a license
also possible, fromthe NumericalRecipesSoftware Web Site (http://www.nr.com).
UnlicensedtransferofNumericalRecipesprogramstoanyotherformat,ortoanycomputer
except one that is specifically licensed, is strictly prohibited. Technical questions,
corrections, and requests for information should be addressed to Numerical Recipes
Software, P.O. Box 243, Cambridge, MA 02238(USA), email “[email protected]”, or fax
781 863-1739.
LibraryofCongressCataloginginPublicationData
NumericalrecipesinC:theartofscientificcomputing/WilliamH.Press
... [et al.]. – 2nd ed.
Includesbibliographicalreferences(p. ) andindex.
ISBN 0-521-43108-5
1. Numericalanalysis–Computerprograms.2. Science–Mathematics–Computerprograms.
3. C(Computerprogramlanguage) I. Press, WilliamH.
QA297.N866 1992
519.4(cid:1)0285(cid:1)53–dc20 92-8876
AcatalogrecordforthisbookisavailablefromtheBritishLibrary.
ISBN 0521431085 Book
ISBN 0521437202 ExamplebookinC
ISBN 0521437245 Cdiskette(IBM3.5(cid:1)(cid:1),1.44M)
ISBN 0521576083 CDROM(IBMPC/Macintosh)
ISBN 0521576075 CDROM(UNIX)
Contents
PrefacetotheSecondEdition xi
PrefacetotheFirstEdition xiv
LicenseInformation xvi
ComputerProgramsbyChapterandSection xix
1 Preliminaries 1
1.0Introduction 1
1.1ProgramOrganizationandControlStructures 5
1.2SomeCConventionsforScientificComputing 15
1.3Error,Accuracy,andStability 28
2 SolutionofLinearAlgebraicEquations 32
2.0Introduction 32
2.1Gauss-JordanElimination 36
2.2GaussianEliminationwithBacksubstitution 41
2.3LUDecompositionandItsApplications 43
2.4TridiagonalandBandDiagonalSystemsofEquations 50
2.5IterativeImprovementofaSolutiontoLinearEquations 55
2.6SingularValueDecomposition 59
2.7SparseLinearSystems 71
2.8VandermondeMatricesandToeplitzMatrices 90
2.9CholeskyDecomposition 96
2.10QRDecomposition 98
2.11IsMatrixInversionanN3Process? 102
3 InterpolationandExtrapolation 105
3.0Introduction 105
3.1PolynomialInterpolationandExtrapolation 108
3.2RationalFunctionInterpolationandExtrapolation 111
3.3CubicSplineInterpolation 113
3.4HowtoSearchanOrderedTable 117
3.5CoefficientsoftheInterpolatingPolynomial 120
3.6InterpolationinTwoorMoreDimensions 123
v
vi Contents
4 IntegrationofFunctions 129
4.0Introduction 129
4.1ClassicalFormulasforEquallySpacedAbscissas 130
4.2ElementaryAlgorithms 136
4.3RombergIntegration 140
4.4ImproperIntegrals 141
4.5GaussianQuadraturesandOrthogonalPolynomials 147
4.6MultidimensionalIntegrals 161
5 EvaluationofFunctions 165
5.0Introduction 165
5.1SeriesandTheirConvergence 165
5.2EvaluationofContinuedFractions 169
5.3PolynomialsandRationalFunctions 173
5.4ComplexArithmetic 176
5.5RecurrenceRelationsandClenshaw’sRecurrenceFormula 178
5.6QuadraticandCubicEquations 183
5.7NumericalDerivatives 186
5.8ChebyshevApproximation 190
5.9DerivativesorIntegralsofaChebyshev-approximatedFunction 195
5.10PolynomialApproximationfromChebyshevCoefficients 197
5.11EconomizationofPowerSeries 198
5.12Pade´ Approximants 200
5.13RationalChebyshevApproximation 204
5.14EvaluationofFunctionsbyPathIntegration 208
6 SpecialFunctions 212
6.0Introduction 212
6.1GammaFunction,BetaFunction,Factorials,BinomialCoefficients 213
6.2IncompleteGammaFunction,ErrorFunction,Chi-Square
ProbabilityFunction,CumulativePoissonFunction 216
6.3ExponentialIntegrals 222
6.4IncompleteBetaFunction,Student’sDistribution,F-Distribution,
CumulativeBinomialDistribution 226
6.5BesselFunctionsofIntegerOrder 230
6.6ModifiedBesselFunctionsofIntegerOrder 236
6.7BesselFunctionsofFractionalOrder,AiryFunctions,Spherical
BesselFunctions 240
6.8SphericalHarmonics 252
6.9FresnelIntegrals,CosineandSineIntegrals 255
6.10Dawson’sIntegral 259
6.11EllipticIntegralsandJacobianEllipticFunctions 261
6.12HypergeometricFunctions 271
7 RandomNumbers 274
7.0Introduction 274
7.1UniformDeviates 275
Contents vii
7.2TransformationMethod: ExponentialandNormalDeviates 287
7.3RejectionMethod: Gamma, Poisson,BinomialDeviates 290
7.4GenerationofRandomBits 296
7.5RandomSequencesBasedonDataEncryption 300
7.6SimpleMonteCarloIntegration 304
7.7Quasi-(thatis,Sub-)RandomSequences 309
7.8AdaptiveandRecursiveMonteCarloMethods 316
8 Sorting 329
8.0Introduction 329
8.1StraightInsertionandShell’sMethod 330
8.2Quicksort 332
8.3Heapsort 336
8.4IndexingandRanking 338
8.5SelectingtheMthLargest 341
8.6DeterminationofEquivalenceClasses 345
9 RootFindingandNonlinearSetsofEquations 347
9.0Introduction 347
9.1BracketingandBisection 350
9.2SecantMethod,FalsePositionMethod,andRidders’Method 354
9.3VanWijngaarden–Dekker–BrentMethod 359
9.4Newton-RaphsonMethodUsingDerivative 362
9.5RootsofPolynomials 369
9.6Newton-RaphsonMethodforNonlinearSystemsofEquations 379
9.7GloballyConvergentMethodsforNonlinearSystemsofEquations 383
10 MinimizationorMaximizationofFunctions 394
10.0Introduction 394
10.1GoldenSectionSearchinOneDimension 397
10.2ParabolicInterpolationandBrent’sMethodinOneDimension 402
10.3One-DimensionalSearchwithFirstDerivatives 405
10.4DownhillSimplexMethodinMultidimensions 408
10.5DirectionSet(Powell’s)MethodsinMultidimensions 412
10.6ConjugateGradientMethodsinMultidimensions 420
10.7VariableMetricMethodsinMultidimensions 425
10.8LinearProgrammingandtheSimplexMethod 430
10.9SimulatedAnnealingMethods 444
11 Eigensystems 456
11.0Introduction 456
11.1JacobiTransformationsofaSymmetricMatrix 463
11.2ReductionofaSymmetricMatrixtoTridiagonalForm:
GivensandHouseholderReductions 469
11.3EigenvaluesandEigenvectorsofaTridiagonalMatrix 475
11.4HermitianMatrices 481
11.5ReductionofaGeneralMatrixtoHessenbergForm 482
viii Contents
11.6TheQRAlgorithmforRealHessenbergMatrices 486
11.7ImprovingEigenvaluesand/orFindingEigenvectorsby
InverseIteration 493
12 FastFourierTransform 496
12.0Introduction 496
12.1FourierTransformofDiscretelySampledData 500
12.2FastFourierTransform(FFT) 504
12.3FFTofRealFunctions,SineandCosineTransforms 510
12.4FFTinTwoorMoreDimensions 521
12.5FourierTransformsofRealDatainTwoandThreeDimensions 525
12.6ExternalStorageorMemory-LocalFFTs 532
13 FourierandSpectralApplications 537
13.0Introduction 537
13.1ConvolutionandDeconvolutionUsingtheFFT 538
13.2CorrelationandAutocorrelationUsingtheFFT 545
13.3Optimal(Wiener)FilteringwiththeFFT 547
13.4PowerSpectrumEstimationUsingtheFFT 549
13.5DigitalFilteringintheTimeDomain 558
13.6LinearPredictionandLinearPredictiveCoding 564
13.7PowerSpectrumEstimationbytheMaximumEntropy
(AllPoles)Method 572
13.8SpectralAnalysisofUnevenlySampledData 575
13.9ComputingFourierIntegralsUsingtheFFT 584
13.10WaveletTransforms 591
13.11NumericalUseoftheSamplingTheorem 606
14 StatisticalDescriptionofData 609
14.0Introduction 609
14.1MomentsofaDistribution:Mean,Variance,Skewness,
andSoForth 610
14.2DoTwoDistributionsHavetheSameMeansorVariances? 615
14.3AreTwoDistributionsDifferent? 620
14.4ContingencyTableAnalysisofTwoDistributions 628
14.5LinearCorrelation 636
14.6NonparametricorRankCorrelation 639
14.7DoTwo-DimensionalDistributionsDiffer? 645
14.8Savitzky-GolaySmoothingFilters 650
15 ModelingofData 656
15.0Introduction 656
15.1LeastSquaresasaMaximumLikelihoodEstimator 657
15.2FittingDatatoaStraightLine 661
15.3Straight-LineDatawithErrorsinBothCoordinates 666
15.4GeneralLinearLeastSquares 671
15.5NonlinearModels 681
Contents ix
15.6ConfidenceLimitsonEstimatedModelParameters 689
15.7RobustEstimation 699
16 IntegrationofOrdinaryDifferentialEquations 707
16.0Introduction 707
16.1Runge-KuttaMethod 710
16.2AdaptiveStepsizeControlforRunge-Kutta 714
16.3ModifiedMidpointMethod 722
16.4RichardsonExtrapolationandtheBulirsch-StoerMethod 724
16.5Second-OrderConservativeEquations 732
16.6StiffSetsofEquations 734
16.7Multistep,Multivalue,andPredictor-CorrectorMethods 747
17 TwoPointBoundaryValueProblems 753
17.0Introduction 753
17.1TheShootingMethod 757
17.2ShootingtoaFittingPoint 760
17.3RelaxationMethods 762
17.4AWorkedExample: SpheroidalHarmonics 772
17.5AutomatedAllocationofMeshPoints 783
17.6HandlingInternalBoundaryConditionsorSingularPoints 784
18 IntegralEquationsandInverseTheory 788
18.0Introduction 788
18.1FredholmEquationsoftheSecondKind 791
18.2VolterraEquations 794
18.3IntegralEquationswithSingularKernels 797
18.4InverseProblemsandtheUseofAPrioriInformation 804
18.5LinearRegularizationMethods 808
18.6Backus-GilbertMethod 815
18.7MaximumEntropyImageRestoration 818
19 PartialDifferentialEquations 827
19.0Introduction 827
19.1Flux-ConservativeInitialValueProblems 834
19.2DiffusiveInitialValueProblems 847
19.3InitialValueProblemsinMultidimensions 853
19.4FourierandCyclicReductionMethodsforBoundary
ValueProblems 857
19.5RelaxationMethodsforBoundaryValueProblems 863
19.6MultigridMethodsforBoundaryValueProblems 871
20 Less-NumericalAlgorithms 889
20.0Introduction 889
20.1DiagnosingMachineParameters 889
20.2GrayCodes 894
x Contents
20.3CyclicRedundancyandOtherChecksums 896
20.4HuffmanCodingandCompressionofData 903
20.5ArithmeticCoding 910
20.6ArithmeticatArbitraryPrecision 915
References 926
AppendixA:TableofPrototypeDeclarations 930
AppendixB:UtilityRoutines 940
AppendixC:ComplexArithmetic 948
IndexofProgramsandDependencies 951
GeneralIndex 965
Preface to the Second Edition
OuraiminwritingtheoriginaleditionofNumericalRecipeswas toprovidea
book that combined general discussion, analytical mathematics, algorithmics, and
actualworkingprograms. Thesuccessofthefirsteditionputsusnowinadifficult,
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.
Itisamixedblessingthatweknowmorenowthanwedidsixyearsago. Then,
weweremakingeducatedguesses,basedonexistingliteratureandourownresearch,
aboutwhichnumericaltechniqueswerethemostimportantandrobust. Now,wehave
thebenefitofdirectfeedbackfromalargereadercommunity. Letterstoouralter-ego
enterprise,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 a
particular field of science or engineering. We value such letters, and digest them
carefully,especiallywhentheypointustospecificreferences intheliterature.
The inevitable result of this input is that this Second Edition of Numerical
Recipes issubstantiallylargerthanitspredecessor, infactabout50%largerbothin
wordsandnumberofincludedprograms(thelatternownumberingwellover300).
“Don’tletthebookgrowinsize,” istheadvice thatwe receivedfromseveral wise
colleagues. We have triedto followthe intended spiritof that advice, even as we
violatetheletterofit. Wehavenotlengthened,orincreasedindifficulty,thebook’s
principaldiscussionsof mainstream topics. Many new topicsare presented at this
same accessible level. Some topics, both from the earlier edition and new to this
one, are nowset insmaller type thatlabelsthem as being“advanced.” The reader
whoignoressuchadvancedsectionscompletelywillnot,wethink,findanylackof
continuity in the shorter volume that results.
Hereare somehighlightsofthenewmaterialinthisSecondEdition:
• a new chapter onintegralequationsand inverse methods
• a detailed treatment of multigrid methods for solving elliptic partial
differential equations
• routines for band diagonal linear systems
• improvedroutinesfor linearalgebra onsparse matrices
• Cholesky and QR decomposition
• orthogonal polynomials and Gaussian quadratures for arbitrary weight
functions
• methods for calculating numerical derivatives
• Pade´ approximants,and rationalChebyshevapproximation
• Bessel functions,and modified Bessel functions,of fractionalorder; 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
• globallyconvergentmethodsforsetsofnonlinearequations
xi
xii PrefacetotheSecondEdition
• simulatedannealingminimizationforcontinuouscontrolspaces
• fastFouriertransform(FFT)forrealdataintwoandthreedimensions
• fast Fourier transform (FFT) usingexternal 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
• fittingstraightlinedata witherrors inbothcoordinates
• a two-dimensionalKolmogorov-Smirnofftest
• the statistical bootstrap method
• embeddedRunge-Kutta-Fehlbergmethodsfordifferentialequations
• high-ordermethods for stiff differential equations
• a new chapter on “less-numerical” algorithms, including Huffman and
arithmeticcoding,arbitraryprecisionarithmetic,andseveralothertopics.
Consultthe Preface tothe FirstEdition, following,or the Table of Contents, for a
list 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;wearegratefulforthese. Inthetext,weattempttogivespecific
attributionforideasthatappear tobe original,andnotknownintheliterature. 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,Sallie
Baliunas, 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, Bradley
Keister,andWilliamGould. WehavebeenhelpedbyNancyLee Snyder’smastery
of a complicated TEX manuscript. We express appreciation to our editors Lauren
Cowles and Alan Harvey at Cambridge University Press, and to our production
editorRussellHahn. Weremain,ofcourse,gratefultotheindividualsacknowledged
in the Preface to the First Edition.
Special acknowledgment is due to programming consultant Seth Finkelstein,
who wrote, rewrote, or influenced many of the routinesin thisbook, as wellas in
its FORTRAN-language twin and the companion Example books. Our project has
benefitedenormouslyfromSeth’stalentfordetecting,andfollowingthetrailof,even
very slightanomalies (often compilerbugs, butoccasionally our errors), and from
hisgoodprogrammingsense. To the extentthatthiseditionofNumericalRecipes
inChasamoregracefuland“C-like”programmingstylethanitspredecessor,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 run-
ning the UNIX operating system, and on a 486/33 PC compatible running
MS-DOS5.0/Windows3.0. (See §1.0 for a list of additional computers used in
