xgmainc.f

Go to the documentation of this file.

      subroutine xgammainc (a, x, result)

c -- jwe, based on DGAMIT.
c
c -- Do a better job than dgami for large values of x.

      double precision a, x, result
      intrinsic exp, log, sqrt, sign, aint
      external dgami, dlngam, d9lgit, d9lgic, d9gmit

C     external dgamr
C     DOUBLE PRECISION DGAMR

      DOUBLE PRECISION aeps, ainta, algap1, alneps, alng, alx,
     $     bot, h, sga, sgngam, sqeps, t, d1mach, d9gmit,
     $     d9lgic, d9lgit, dlngam, dgami 

      LOGICAL first

      SAVE alneps, sqeps, bot, first

      DATA first /.true./

      if (x .eq. 0.0d0) then

        if (a .eq. 0.0d0) then
          result = 1.0d0
        else
          result = 0.0d0
        endif

      else

      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', 'XGMAINC', '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
C      IF (X.GT.0.D0) GO TO 20
C      DGAMIT = 0.0D0
C      IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
C      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)
C      DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
      result = exp(a*alx + log(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
C      DGAMIT = EXP (T)
      result = exp(a*alx + 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', 'XGMAINC', 'RESULT LT HALF PRECISION', 1,
     +   1)
C
C 50   T = -A*ALX + LOG(ABS(H))
C      IF (T.LT.BOT) CALL XERCLR
C      DGAMIT = SIGN (EXP(T), H)
 50   result = h
      RETURN
C
C 60   T = T - A*ALX
 60   IF (t.LT.bot) CALL xerclr
      result = -sga * sgngam * exp(t)
      RETURN

      endif
      return
      end