dgamit.f

Go to the documentation of this file.

*DECK DGAMIT
      DOUBLE PRECISION FUNCTION dgamit (A, X)
C***BEGIN PROLOGUE  DGAMIT
C***PURPOSE  Calculate Tricomi's form of the incomplete Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      DOUBLE PRECISION (GAMIT-S, DGAMIT-D)
C***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
C             SPECIAL FUNCTIONS, TRICOMI
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C   Evaluate Tricomi's incomplete Gamma function defined by
C
C   DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
C              T**(A-1.)
C
C   for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
C   GAMMA(X) is the complete gamma function of X.
C
C   DGAMIT is evaluated for arbitrary real values of A and for non-
C   negative values of X (even though DGAMIT is defined for X .LT.
C   0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite,
C   which is a fatal error.
C
C   The function and both arguments are DOUBLE PRECISION.
C
C   A slight deterioration of 2 or 3 digits accuracy will occur when
C   DGAMIT is very large or very small in absolute value, because log-
C   arithmic variables are used.  Also, if the parameter  A  is very
C   close to a negative integer (but not a negative integer), there is
C   a loss of accuracy, which is reported if the result is less than
C   half machine precision.
C
C***REFERENCES  W. Gautschi, A computational procedure for incomplete
C                 gamma functions, ACM Transactions on Mathematical
C                 Software 5, 4 (December 1979), pp. 466-481.
C               W. Gautschi, Incomplete gamma functions, Algorithm 542,
C                 ACM Transactions on Mathematical Software 5, 4
C                 (December 1979), pp. 482-489.
C***ROUTINES CALLED  D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS,
C                    DLNGAM, XERCLR, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
C***END PROLOGUE  DGAMIT
      DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
     1  bot, h, sga, sgngam, sqeps, t, d1mach, dgamr, d9gmit, d9lgit,
     2  dlngam, d9lgic
      LOGICAL FIRST
      SAVE alneps, sqeps, bot, first
      DATA first /.true./
C***FIRST EXECUTABLE STATEMENT  DGAMIT
      IF (first) THEN
         alneps = -log(d1mach(3))
         sqeps = sqrt(d1mach(4))
         bot = log(d1mach(1))
      ENDIF
      first = .false.
C
      IF (x .LT. 0.d0) CALL xermsg('SLATEC', 'DGAMIT', 'X IS NEGATIVE'
     +   , 2, 2)
C
      IF (x.NE.0.d0) alx = log(x)
      sga = 1.0d0
      IF (a.NE.0.d0) sga = sign(1.0d0, a)
      ainta = aint(a + 0.5d0*sga)
      aeps = a - ainta
C
      IF (x.GT.0.d0) go to 20
      dgamit = 0.0d0
      IF (ainta.GT.0.d0 .OR. aeps.NE.0.d0) dgamit = dgamr(a+1.0d0)
      RETURN
C
 20   IF (x.GT.1.d0) go to 30
      IF (a.GE.(-0.5d0) .OR. aeps.NE.0.d0) CALL dlgams(a+1.0d0, algap1,
     1  sgngam)
      dgamit = d9gmit(a, x, algap1, sgngam, alx)
      RETURN
C
 30   IF (a.LT.x) go to 40
      t = d9lgit(a, x, dlngam(a+1.0d0))
      IF (t.LT.bot) CALL xerclr
      dgamit = exp(t)
      RETURN
C
 40   alng = d9lgic(a, x, alx)
C
C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
C
      h = 1.0d0
      IF (aeps.EQ.0.d0 .AND. ainta.LE.0.d0) go to 50
C
      CALL dlgams(a+1.0d0, algap1, sgngam)
      t = log(abs(a)) + alng - algap1
      IF (t.GT.alneps) go to 60
C
      IF (t.GT.(-alneps)) h = 1.0d0 - sga * sgngam * exp(t)
      IF (abs(h).GT.sqeps) go to 50
C
      CALL xerclr
      CALL xermsg('SLATEC', 'DGAMIT', 'RESULT LT HALF PRECISION', 1,
     +   1)
C
 50   t = -a*alx + log(abs(h))
      IF (t.LT.bot) CALL xerclr
      dgamit = sign(exp(t), h)
      RETURN
C
 60   t = t - a*alx
      IF (t.LT.bot) CALL xerclr
      dgamit = -sga * sgngam * exp(t)
      RETURN
C
      END