xsgmainc.f

Go to the documentation of this file.

      subroutine xsgammainc (a, x, result)

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

      real a, x, result
      intrinsic exp, log, sqrt, sign, aint
      external gami, alngam, r9lgit, r9lgic, r9gmit

C     external gamr
C     real GAMR

      REAL AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
     $     bot, h, sga, sgngam, sqeps, t, r1mach, r9gmit,
     $     r9lgic, r9lgit, alngam, gami

      LOGICAL FIRST

      SAVE alneps, sqeps, bot, first

      DATA first /.true./

      if (x .eq. 0.0e0) then

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

      else

      IF (first) THEN
         alneps = -log(r1mach(3))
         sqeps = sqrt(r1mach(4))
         bot = log(r1mach(1))
      ENDIF
      first = .false.
C
      IF (x .LT. 0.e0) CALL xermsg('SLATEC', 'XGMAINC', 'X IS NEGATIVE'
     +   , 2, 2)
C
      IF (x.NE.0.e0) alx = log(x)
      sga = 1.0e0
      IF (a.NE.0.e0) sga = sign(1.0e0, a)
      ainta = aint(a + 0.5e0*sga)
      aeps = a - ainta
C
C      IF (X.GT.0.E0) GO TO 20
C      GAMIT = 0.0E0
C      IF (AINTA.GT.0.E0 .OR. AEPS.NE.0.E0) GAMIT = GAMR(A+1.0E0)
C      RETURN
C
 20   IF (x.GT.1.e0) go to 30
      IF (a.GE.(-0.5e0) .OR. aeps.NE.0.e0) CALL algams(a+1.0e0, algap1,
     1  sgngam)
C      GAMIT = R9GMIT (A, X, ALGAP1, SGNGAM, ALX)
      result = exp(a*alx + log(r9gmit(a, x, algap1, sgngam, alx)))
      RETURN
C
 30   IF (a.LT.x) go to 40
      t = r9lgit(a, x, alngam(a+1.0e0))
      IF (t.LT.bot) CALL xerclr
C      GAMIT = EXP (T)
      result = exp(a*alx + t)
      RETURN
C
 40   alng = r9lgic(a, x, alx)
C
C EVALUATE GAMIT IN TERMS OF LOG (DGAMIC (A, X))
C
      h = 1.0e0
      IF (aeps.EQ.0.e0 .AND. ainta.LE.0.e0) go to 50
C
      CALL algams(a+1.0e0, algap1, sgngam)
      t = log(abs(a)) + alng - algap1
      IF (t.GT.alneps) go to 60
C
      IF (t.GT.(-alneps)) h = 1.0e0 - 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      GAMIT = 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