function r9gmic (a, x, alx) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. c c compute the complementary incomplete gamma function for a near c a negative integer and for small x. c external alngam, r1mach data euler / .5772156649 015329 e0 / data eps, bot / 2*0.0 / c if (eps.eq.0.0) eps = 0.5*r1mach(3) if (bot.eq.0.0) bot = alog(r1mach(1)) c if (a.gt.0.0) call seteru ( 1 41hr9gmic a must be near a negative integer, 41, 2, 2) if (x.le.0.0) call seteru (25hr9gmic x must be gt zero, 25, 3, 2) c ma = a - 0.5 fm = -ma m = -ma c te = 1.0 t = 1.0 s = t do 20 k=1,200 fkp1 = k + 1 te = -x*te/(fm+fkp1) t = te/fkp1 s = s + t if (abs(t).lt.eps*s) go to 30 20 continue call seteru (57hr9gmic no convergence in 200 terms of continued f 1raction, 57, 4, 2) c 30 r9gmic = -alx - euler + x*s/(fm+1.0) if (m.eq.0) return c if (m.eq.1) r9gmic = -r9gmic - 1.0 + 1.0/x if (m.eq.1) return c te = fm t = 1.0 s = t mm1 = m - 1 do 40 k=1,mm1 fk = k te = -x*te/fk t = te/(fm-fk) s = s + t if (abs(t).lt.eps*abs(s)) go to 50 40 continue c 50 do 60 k=1,m r9gmic = r9gmic + 1.0/float(k) 60 continue c sgng = 1.0 if (mod(m,2).eq.1) sgng = -1.0 alng = alog(r9gmic) - alngam(fm+1.0) c r9gmic = 0.0 if (alng.gt.bot) r9gmic = sgng*exp(alng) if (s.ne.0.0) r9gmic = r9gmic + sign (exp(-fm*alx+alog(abs(s)/fm)) 1 , s) c if (r9gmic.eq.0.0 .and. s.eq.0.0) call seteru ( 1 25hr9gmic result underflows, 25, 1, 0) return c end