f6-1

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 6.1Gamma, Beta,andRelatedFunctions 207Sample 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).same parameters) the formula (6.1.5)and bound on /epsilon1apply for the complexgamma function, everywhere in the half complex plane Re z> 0. It is better to implement ln Γ(x)than Γ(x), since the latter will overflowmany computers’ floating-point representation at quite modest values of x. Often the gammafunctionisusedincalculationswherethelargevaluesof Γ(x)aredividedby otherlargenumbers,withtheresultbeingaperfectlyordinaryvalue. Suchoperations wouldnormallybecodedas subtractionoflogarithms. With (6.1.5)in hand,wecan compute the logarithm of the gamma function with two calls to a logarithm and 25 orsoarithmeticoperations. Thismakesitnotmuchmoredifficultthanotherbuilt-in functions that we take for granted, such as sinxorex: FUNCTION gammln(xx) REAL gammln,xx Returns the value ln[Γ( xx)]forxx >0. INTEGER j DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) Internal arithmetic willbedone indouble precision, anicety that youcanomitiffive-figureaccuracy is good enough. SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, * 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,* -.5395239384953d-5,2.5066282746310005d0/ x=xx y=xtmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0do 11j=1,6 y=y+1.d0 ser=ser+cof(j)/y enddo 11 gammln=tmp+log(stp*ser/x) return END How shall we write a routine for the factorial function n!? Generally the factorial function will be called for small integer values (for large values it will overflowanyway!),andinmostapplicationsthesameintegervaluewillbecalledfor manytimes. Itis aprofligatewaste ofcomputertimetocall exp(gammln(n+1.0)) for each required factorial. Better to go back to basics, holding gammlnin reserve for unlikely calls: FUNCTION factrl(n) INTEGER nREAL factrl C USES gammln Returns the value n!as a floating-point number. INTEGER j,ntopREAL a(33),gammln Table to be filled in only as required. SAVE ntop,a DATA ntop,a(1)/0,1./ Table initialized with 0!only. if (n.lt.0) then pause ’negative factorial in factrl’ else if (n.le.ntop) then Already in table. factrl=a(n+1) else if (n.le.32) then Fill in table up to desired value. do 11j=ntop+1,n 208 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).a(j+1)=j*a(j) enddo 11 ntop=nfactrl=a(n+1) else Larger value than size of table is required. Actually, this big a value is going to overflow on many computers, but no harm in trying.factrl=exp(gammln(n+1.)) endifreturnEND A useful point is that factrlwill beexactfor the smaller values of n, since floating-pointmultipliesonsmallintegersareexactonallcomputers. Thisexactness will not hold if we turn to the logarithmof the factorials. For binomial coefficients, however, we must do exactly this, since the individual factorials in a binomial coefficient will overflow long before the coefficient itself will. The binomial coefficient is defined by /parenleftbiggn k/parenrightbigg =n! k!(n−k)!0≤k≤n (6.1.6 ) FUNCTION bico(n,k) INTEGER k,nREAL bico C USES factln Returns the binomial coefficient/parenleftbign k/parenrightbigas a floating-point number. REAL factlnbico=nint(exp(factln(n)-factln(k)-factln(n-k)))returnThe nearest-integer function cleans up roundoff error for smaller values of nandk. END which uses FUNCTION factln(n) INTEGER nREAL factln C USES gammln Returns ln(n!). REAL a(100),gammlnSAVE a DATA a/100*-1./ Initializethe tableto negative values. if (n.lt.0) pause ’negative factorial in factln’if (n.le.99) then In range of the table. if (a(n+1).lt.0.) a(n+1)=gammln(n+1.) Ifnot already inthe table, put it in. factln=a(n+1) else factln=gammln(n+1.) Out of range of the table. endif returnEND 6.2IncompleteGammaFunction 209Sample 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).If your problem requires a series of related binomial coefficients, a good idea is to use recurrence relations, for example /parenleftbiggn+1 k/parenrightbigg =n+1 n−k+1/parenleftbiggn k/parenrightbigg =/parenleftbiggn k/parenrightbigg +/parenleftbiggn k−1/parenrightbigg /parenleftbiggn k+1/parenrightbigg =n−k k+1/parenleftbiggn k/parenrightbigg (6.1.7 ) Finally, turning away from the combinatorial functions with integer valued arguments, we come to the beta function, B(z,w)=B(w, z)=/integraldisplay1 0tz−1(1−t)w−1dt (6.1.8 ) which is related to the gamma function by B(z,w)=Γ(z)Γ(w) Γ(z+w)(6.1.9 ) hence FUNCTION beta(z,w) REAL beta,w,z C USES gammln Returns the value of the beta function B(z, w ). REAL gammlnbeta=exp(gammln(z)+gammln(w)-gammln(z+w))return END 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), Chapter 6. Lanczos, C. 1964, SIAM Journal on Numerical Analysis , ser. B, vol. 1, pp. 86–96. [1] 6.2 Incomplete Gamma Function, Error Function, Chi-Square Probability Function, Cumulative Poisson Function The incomplete gamma function is defined by P(a, x)≡γ(a, x) Γ(a)≡1 Γ(a)/integraldisplayx 0e−tta−1dt (a> 0) ( 6.2.1 )