function gamic (a, x) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. c c evaluate the complementary incomplete gamma function c c gamic = integral from t = x to infinity of exp(-t) * t**(a-1.) . c c gamic is evaluated for arbitrary real values of a and for non-negative c values of x (even though gamic is defined for x .lt. 0.0), except that c for x = 0 and a .le. 0.0, gamic is undefined. c c a slight deterioration of 2 or 3 digits accuracy will occur when c gamic is very large or very small in absolute value, because log- c arithmic variables are used. also, if the parameter a is very close c to a negative integer (but not a negative integer), there is a loss c of accuracy, which is reported if the result is less than half c machine precision. c c ref. -- w. gautschi, an evaluation procedure for incomplete gamma c functions, acm trans. math. software. c external alngam, r1mach, r9gmic, r9gmit, r9lgic, 1 r9lgit data eps, sqeps, alneps, bot / 4*0.0 / c if (eps.ne.0.) go to 10 eps = 0.5*r1mach(3) sqeps = sqrt(r1mach(4)) alneps = -alog(r1mach(3)) bot = alog(r1mach(1)) c 10 if (x.lt.0.0) call seteru (21hgamic x is negative, 21, 2, 2) c if (x.gt.0.0) go to 20 if (a.le.0.0) call seteru ( 1 46hgamic x = 0 and a le 0 so gamic is undefined, 46, 3, 2) c gamic = exp (alngam(a+1.0) - alog(a)) return c 20 alx = alog(x) sga = 1.0 if (a.ne.0.0) sga = sign (1.0, a) ma = a + 0.5*sga aeps = a - float(ma) c izero = 0 if (x.ge.1.0) go to 60 c if (a.gt.0.5 .or. abs(aeps).gt.0.001) go to 50 fm = -ma e = 2.0 if (fm.gt.1.0) e = 2.0*(fm+2.0)/(fm*fm-1.0) e = e - alx*x**(-0.001) if (e*abs(aeps).gt.eps) go to 50 c gamic = r9gmic (a, x, alx) return c 50 call algams (a+1.0, algap1, sgngam) gstar = r9gmit (a, x, algap1, sgngam, alx) if (gstar.eq.0.0) izero = 1 if (gstar.ne.0.0) alngs = alog (abs(gstar)) if (gstar.ne.0.0) sgngs = sign (1.0, gstar) go to 70 c 60 if (a.lt.x) gamic = exp (r9lgic(a, x, alx)) if (a.lt.x) return c sgngam = 1.0 algap1 = alngam (a+1.0) sgngs = 1.0 alngs = r9lgit (a, x, algap1) c c evaluation of gamic(a,x) in terms of tricomi-s incomplete gamma fn. c 70 h = 1.0 if (izero.eq.1) go to 80 c t = a*alx + alngs if (t.gt.alneps) go to 90 if (t.gt.(-alneps)) h = 1.0 - sgngs*exp(t) c if (abs(h).lt.sqeps) call erroff if (abs(h).lt.sqeps) call seteru ( 1 32hgamic result lt half precision, 32, 1, 1) c 80 sgng = sign (1.0, h) * sga * sgngam t = alog(abs(h)) + algap1 - alog(abs(a)) if (t.lt.bot) call erroff gamic = sgng * exp(t) return c 90 sgng = -sgngs * sga * sgngam t = t + algap1 - alog(abs(a)) if (t.lt.bot) call erroff gamic = sgng * exp(t) return c end