Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 56
Текст из файла (страница 56)
Actually, this big a value is going to overflowon many computers, but no harm in trying.while (ntop<n) {Fill in table up to desired value.j=ntop++;a[ntop]=a[j]*ntop;}return a[n];}6.1 Gamma, Beta, and Related Functions215A useful point is that factrl will be exact for the smaller values of n, sincefloating-point multiplies on small integers are exact on all computers. This exactnesswill not hold if we turn to the logarithm of the factorials.
For binomial coefficients,however, we must do exactly this, since the individual factorials in a binomialcoefficient will overflow long before the coefficient itself will.The binomial coefficient is defined by n!n=k!(n − k)!k0≤k≤n(6.1.6)#include <math.h>float bico(int n, int k), Returns the binomial coefficient nas a floating-point number.k{float factln(int n);return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));The floor function cleans up roundoff error for smaller values of n and k.}which usesfloat factln(int n)Returns ln(n!).{float gammln(float xx);void nrerror(char error_text[]);static float a[101];A static array is automatically initialized to zero.if (n < 0) nrerror("Negative factorial in routine factln");if (n <= 1) return 0.0;if (n <= 100) return a[n] ? a[n] : (a[n]=gammln(n+1.0)); In range of table.else return gammln(n+1.0);Out of range of table.}If your problem requires a series of related binomial coefficients, a good ideais to use recurrence relations, for example n+1n+1nnn==+kn−k+1 kkk−1 n−k nn=k+1 kk+1(6.1.7)Finally, turning away from the combinatorial functions with integer valuedarguments, we come to the beta function,&1B(z, w) = B(w, z) =0tz−1 (1 − t)w−1 dt(6.1.8)216Chapter 6.Special Functionswhich is related to the gamma function byB(z, w) =Γ(z)Γ(w)Γ(z + w)(6.1.9)hence#include <math.h>float beta(float z, float w)Returns the value of the beta function B(z, w).{float gammln(float xx);return exp(gammln(z)+gammln(w)-gammln(z+w));}CITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A.
1964, Handbook of Mathematical Functions, Applied Mathematics 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, ErrorFunction, Chi-Square Probability Function,Cumulative Poisson FunctionThe incomplete gamma function is defined byP (a, x) ≡1γ(a, x)≡Γ(a)Γ(a)&xe−t ta−1 dt(a > 0)(6.2.1)0It has the limiting valuesP (a, 0) = 0P (a, ∞) = 1and(6.2.2)The incomplete gamma function P (a, x) is monotonic and (for a greater than one orso) rises from “near-zero”to “near-unity” in a range of x centered on about a − 1,√and of width about a (see Figure 6.2.1).The complement of P (a, x) is also confusingly called an incomplete gammafunction,Q(a, x) ≡ 1 − P (a, x) ≡1Γ(a, x)≡Γ(a)Γ(a)&∞xe−t ta−1 dt(a > 0) (6.2.3)2176.2 Incomplete Gamma Function1.0incomplete gamma function P(a,x)0.51.0.8a = 3.0.6a = 10.4.2002468101214xFigure 6.2.1.The incomplete gamma function P (a, x) for four values of a.It has the limiting valuesQ(a, 0) = 1andQ(a, ∞) = 0(6.2.4)The notations P (a, x), γ(a, x), and Γ(a, x) are standard; the notation Q(a, x) isspecific to this book.There is a series development for γ(a, x) as follows:γ(a, x) = e−x xa∞Γ(a)xnΓ(a+1+n)n=0(6.2.5)One does not actually need to compute a new Γ(a + 1 + n) for each n; one ratheruses equation (6.1.3) and the previous coefficient.A continued fraction development for Γ(a, x) is1 1−a 1 2−a 2···(x > 0)(6.2.6)Γ(a, x) = e−x xax+ 1+ x+ 1+ x+It is computationally better to use the even part of (6.2.6), which converges twiceas fast (see §5.2):1 · (1 − a)2 · (2 − a)1−x a···(x > 0)Γ(a, x) = e xx+1−a− x+3−a− x+5−a−(6.2.7)It turns out that (6.2.5) converges rapidly for x less than about a + 1, while(6.2.6) or (6.2.7) converges rapidly for x greater than about a + 1.
In these respective218Chapter 6.Special Functions√regimes each requires at most a few times a terms to converge, and this manyonly near x = a, where the incomplete gamma functions are varying most rapidly.Thus (6.2.5) and (6.2.7) together allow evaluation of the function for all positivea and x. An extra dividend is that we never need compute a function value nearzero by subtracting two nearly equal numbers.
The higher-level functions that returnP (a, x) and Q(a, x) arefloat gammp(float a, float x)Returns the incomplete gamma function P (a, x).{void gcf(float *gammcf, float a, float x, float *gln);void gser(float *gamser, float a, float x, float *gln);void nrerror(char error_text[]);float gamser,gammcf,gln;if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp");if (x < (a+1.0)) {Use the series representation.gser(&gamser,a,x,&gln);return gamser;} else {Use the continued fraction representationgcf(&gammcf,a,x,&gln);return 1.0-gammcf;and take its complement.}}float gammq(float a, float x)Returns the incomplete gamma function Q(a, x) ≡ 1 − P (a, x).{void gcf(float *gammcf, float a, float x, float *gln);void gser(float *gamser, float a, float x, float *gln);void nrerror(char error_text[]);float gamser,gammcf,gln;if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq");if (x < (a+1.0)) {Use the series representationgser(&gamser,a,x,&gln);return 1.0-gamser;and take its complement.} else {Use the continued fraction representation.gcf(&gammcf,a,x,&gln);return gammcf;}}The argument gln is set by both the series and continued fraction proceduresto the value ln Γ(a); the reason for this is so that it is available to you if you want tomodify the above two procedures to give γ(a, x) and Γ(a, x), in addition to P (a, x)and Q(a, x) (cf.
equations 6.2.1 and 6.2.3).The functions gser and gcf which implement (6.2.5) and (6.2.7) are#include <math.h>#define ITMAX 100#define EPS 3.0e-7void gser(float *gamser, float a, float x, float *gln)Returns the incomplete gamma function P (a, x) evaluated by its series representation as gamser.Also returns ln Γ(a) as gln.{float gammln(float xx);6.2 Incomplete Gamma Function219void nrerror(char error_text[]);int n;float sum,del,ap;*gln=gammln(a);if (x <= 0.0) {if (x < 0.0) nrerror("x less than 0 in routine gser");*gamser=0.0;return;} else {ap=a;del=sum=1.0/a;for (n=1;n<=ITMAX;n++) {++ap;del *= x/ap;sum += del;if (fabs(del) < fabs(sum)*EPS) {*gamser=sum*exp(-x+a*log(x)-(*gln));return;}}nrerror("a too large, ITMAX too small in routine gser");return;}}#include <math.h>#define ITMAX 100#define EPS 3.0e-7#define FPMIN 1.0e-30Maximum allowed number of iterations.Relative accuracy.Number near the smallest representablefloating-point number.void gcf(float *gammcf, float a, float x, float *gln)Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation as gammcf.
Also returns ln Γ(a) as gln.{float gammln(float xx);void nrerror(char error_text[]);int i;float an,b,c,d,del,h;*gln=gammln(a);b=x+1.0-a;Set up for evaluating continued fractionby modified Lentz’s method (§5.2)c=1.0/FPMIN;with b0 = 0.d=1.0/b;h=d;for (i=1;i<=ITMAX;i++) {Iterate to convergence.an = -i*(i-a);b += 2.0;d=an*d+b;if (fabs(d) < FPMIN) d=FPMIN;c=b+an/c;if (fabs(c) < FPMIN) c=FPMIN;d=1.0/d;del=d*c;h *= del;if (fabs(del-1.0) < EPS) break;}if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");*gammcf=exp(-x+a*log(x)-(*gln))*h;Put factors in front.}220Chapter 6.Special FunctionsError FunctionThe error function and complementary error function are special cases of theincomplete gamma function, and are obtained moderately efficiently by the aboveprocedures. Their definitions are2erf(x) = √π&x2e−t dt(6.2.8)0and2erfc(x) ≡ 1 − erf(x) = √π&∞2e−t dt(6.2.9)xThe functions have the following limiting values and symmetries:erf(0) = 0erfc(0) = 1erf(∞) = 1erfc(∞) = 0erf(−x) = −erf(x)erfc(−x) = 2 − erfc(x)(6.2.10)(6.2.11)They are related to the incomplete gamma functions by1 2,xerf(x) = P21erfc(x) = Q , x22and(x ≥ 0)(6.2.12)(x ≥ 0)(6.2.13)We’ll put an extra “f” into our routine names to avoid conflicts with names alreadyin some C libraries:float erff(float x)Returns the error function erf(x).{float gammp(float a, float x);return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);}float erffc(float x)Returns the complementary error function erfc(x).{float gammp(float a, float x);float gammq(float a, float x);return x < 0.0 ? 1.0+gammp(0.5,x*x) : gammq(0.5,x*x);}If you care to do so, you √can easily remedy the minor inefficiency in erff anderffc, namely that Γ(0.5) = π is computed unnecessarily when gammp or gammqis called.
Before you do that, however, you might wish to consider the followingroutine, based on Chebyshev fitting to an inspired guess as to the functional form:2216.2 Incomplete Gamma Function#include <math.h>float erfcc(float x)Returns the complementary error function erfc(x) with fractional error everywhere less than1.2 × 10−7 .{float t,z,ans;z=fabs(x);t=1.0/(1.0+0.5*z);ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))));return x >= 0.0 ? ans : 2.0-ans;}There are also some functions of two variables that are special cases of theincomplete gamma function:Cumulative Poisson Probability FunctionPx (< k), for positive x and integer k ≥ 1, denotes the cumulative Poissonprobability function.
It is defined as the probability that the number of Poissonrandom events occurring will be between 0 and k − 1 inclusive, if the expected meannumber is x. It has the limiting valuesPx (< 1) = e−xPx (< ∞) = 1(6.2.14)Its relation to the incomplete gamma function is simplyPx (< k) = Q(k, x) = gammq (k, x)(6.2.15)Chi-Square Probability FunctionP (χ2 |ν) is defined as the probability that the observed chi-square for a correctmodel should be less than a value χ2 . (We will discuss the use of this function inChapter 15.) Its complement Q(χ2 |ν) is the probability that the observed chi-squarewill exceed the value χ2 by chance even for a correct model.
In both cases ν is aninteger, the number of degrees of freedom. The functions have the limiting valuesP (0|ν) = 0Q(0|ν) = 1P (∞|ν) = 1Q(∞|ν) = 0(6.2.16)(6.2.17)and the following relation to the incomplete gamma functions,ν χ2ν χ2,= gammp,2 22 22ν χ2ν χ,= gammq,Q(χ2 |ν) = Q2 22 2P (χ2 |ν) = P(6.2.18)(6.2.19)222Chapter 6.Special FunctionsCITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapters 6, 7, and 26.Pearson, K. (ed.) 1951, Tables of the Incomplete Gamma Function (Cambridge: CambridgeUniversity Press).6.3 Exponential IntegralsThe standard definition of the exponential integral is&∞En (x) =1e−xtdt,tnx > 0,n = 0, 1, .
. .(6.3.1)The function defined by the principal value of the integral&Ei(x) = −∞−xe−tdt =t&x−∞etdt,tx>0(6.3.2)is also called an exponential integral. Note that Ei(−x) is related to −E1 (x) byanalytic continuation.The function En (x) is a special case of the incomplete gamma functionEn (x) = xn−1 Γ(1 − n, x)(6.3.3)We can therefore use a similar strategy for evaluating it.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
