Prev Next dgamlm.f Headings

      subroutine dgamlm(xmin,xmax)
c***begin prologue  dgamlm
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c7a,r2
c***keywords  complete gamma function,double precision,gamma function,
c             limits,special function
c***author  fullerton, w., (lanl)
c***purpose  computes the d.p. minimum and maximum bounds for x in
c            gamma(x).
c calculate the minimum and maximum legal bounds for x in gamma(x).
c xmin and xmax are not the only bounds, but they are the only non-
c trivial ones to calculate.
c             output arguments --
c xmin   double precision minimum legal value of x in gamma(x).  any
c        smaller value of x might result in underflow.
c xmax   double precision maximum legal value of x in gamma(x).  any
c        larger value of x might cause overflow.
c***references  (none)
c***routines called  d1mach,xerror
c***end prologue  dgamlm
      double precision xmin, xmax, alnbig, alnsml, xln, xold, d1mach
c***first executable statement  dgamlm
      alnsml = dlog(d1mach(1))
      xmin = -alnsml
      do 10 i=1,10
        xold = xmin
        xln = dlog(xmin)
        xmin = xmin - xmin*((xmin+0.5d0)*xln - xmin - 0.2258d0 + alnsml)
     1    / (xmin*xln+0.5d0)
        if (dabs(xmin-xold).lt.0.005d0) go to 20
 10   continue
      call xerror ( 'dgamlm  unable to find xmin', 27, 1, 2)
 20   xmin = -xmin + 0.01d0
      alnbig = dlog (d1mach(2))
      xmax = alnbig
      do 30 i=1,10
        xold = xmax
        xln = dlog(xmax)
        xmax = xmax - xmax*((xmax-0.5d0)*xln - xmax + 0.9189d0 - alnbig)
     1    / (xmax*xln-0.5d0)
        if (dabs(xmax-xold).lt.0.005d0) go to 40
 30   continue
      call xerror ( 'dgamlm  unable to find xmax', 27, 2, 2)
 40   xmax = xmax - 0.01d0
      xmin = dmax1 (xmin, -xmax+1.d0)

Input File: omh/dgamlm.f.omh