f6-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).Chapter 6. Special Functions 6.0 Introduction There is nothing particularly special about a special function , except that some person in authority or textbook writer (not the same thing!) has decided to bestow the moniker. Special functions are sometimes called higher transcendental functions (higherthanwhat?)or functionsofmathematicalphysics (buttheyoccurin other fields also) or functions that satisfy certain frequently occurring second-order differentialequations (butnot all special functionsdo). One mightsimplycall them “useful functions” and let it go at that; it is surely only a matter of taste which functions we have chosen to include in this chapter. Goodcommerciallyavailableprogramlibraries,suchasNAGorIMSL,contain routinesforanumberofspecialfunctions. Theseroutinesareintendedforuserswho will have no idea what goes on inside them. Such state of the art “black boxes” are oftenverymessythings,fullofbranchestocompletelydifferentmethodsdependingon the value of the calling arguments. Black boxes have, or should have, careful control of accuracy, to some stated uniform precision in all regimes. We will not be quite so fastidious in our examples, in part because we want to illustrate techniques from Chapter 5, and in part because we wantyou to understand what goes on in the routines presented. Some of our routines have an accuracy parameter that can be made as small as desired, while others (especially those involving polynomial fits) give only a certain accuracy, one that we believe serviceable (typically six significant figures or more). We do notcertify that the routines are perfect black boxes. We do hope that, if you ever encounter trouble in a routine, you will be able to diagnose and correct the problem on the basis of the information that we have given. In short, the special function routines of this chapter are meant to be used — we use them all the time — but we also want you to be prepared to understand their inner workings. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York) [full of useful numerical approximations to a great varietyof functions]. IMSL Sfun/Library Users Manual (IMSL Inc., 2500 CityWest Boulevard, Houston TX 77042). NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.), Chapter S. 205 206 Chapter6. SpecialFunctionsSample 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).Hart, J.F., et al. 1968, Computer Approximations (New York: Wiley). Hastings,C.1955, ApproximationsforDigitalComputers (Princeton:PrincetonUniversityPress). Luke,Y.L.1975, MathematicalFunctionsandTheirApproximations (NewYork:AcademicPress). 6.1 GammaFunction,BetaFunction,Factorials, Binomial Coefficients The gamma function is defined by the integral Γ(z)=/integraldisplay∞ 0tz−1e−tdt (6.1.1 ) When the argument zis an integer, the gamma functionis just the familiar factorial function, but offset by one, n!=Γ ( n+1 ) ( 6.1.2 ) The gamma function satisfies the recurrence relation Γ(z+1 )= zΓ(z)( 6.1.3 ) Ifthefunctionis knownforarguments z>1or,moregenerally,inthehalfcomplex planeRe (z)>1itcanbeobtainedfor z<1orRe (z)<1bythereflectionformula Γ(1−z)=π Γ(z)s i n ( πz)=πz Γ(1 + z)s i n ( πz)(6.1.4 ) Notice that Γ(z)has a pole at z=0, and at all negativeinteger values of z. There are a variety of methods in use for calculating the function Γ(z) numerically, but none is quite as neat as the approximation derived by Lanczos [1]. This scheme is entirely specific to the gamma function, seemingly plucked from thin air. We will not attempt to derive the approximation, but only state theresultingformula: Forcertainintegerchoicesof γandN,andforcertaincoefficients c 1,c2,...,c N, the gamma function is given by Γ(z+1 )=( z+γ+1 2)z+1 2e−(z+γ+1 2) ×√ 2π/bracketleftbigg c0+c1 z+1+c2 z+2+···+cN z+N+/epsilon1/bracketrightbigg (z>0)(6.1.5 ) You can see that this is a sort of take-off on Stirling’s approximation, but with a series of corrections that take into account the first few poles in the left complex plane. Theconstant c0isverynearlyequalto1. Theerrortermisparametrizedby /epsilon1. Forγ=5,N=6,andacertainsetof c’s, theerroris smallerthan |/epsilon1|<2×10−10. Impressed? If not, then perhaps you will be impressed by the fact that (with these