      SUBROUTINE GAMINV(A,X,X0,P,Q,IERR)
C     ******************************************************************
C     INVERSE INCOMPLETE GAMMA RATIO FUNCTION
C     ******************************************************************
      REAL P0(4)
C     -------------------
      DATA A0/3.31125922108741/, A1/11.6616720288968/,
     *     A2/4.28342155967104/, A3/.213623493715853/
      DATA B1/6.61053765625462/, B2/6.40691597760039/,
     *     B3/1.27364489782223/, B4/.036117081018842/
      DATA C /.577215664901533/
C     -------------------
      DATA TOL/5.E-5/
      DATA P0(1)/.15/, P0(2)/.05/, P0(3)/.01/, P0(4)/.001/
C     -------------------
C     ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS.
C            E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0.
C            XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE
C            LARGEST POSITIVE NUMBER.
C
                   E = SPMPAR(1)
                   XMIN = SPMPAR(2)
                   XMAX = SPMPAR(3)
C     -------------------
      X = 0.0
      IF (A .LE. 0.0) GO TO 500
      T = DBLE(P) + DBLE(Q) - 1.D0
      IF (ABS(T) .GT. E) GO TO 510
      IF (X0 .LT. 0.0) GO TO 520
C
      IERR = 0
      IF (P .EQ. 0.0) RETURN
      IF (Q .EQ. 0.0) GO TO 400
      IF (A .EQ. 1.0) GO TO 410
C
      E2 = 2.0*E
      AMAX = 0.5E-9/(E*E)
      XN = X0
      IF (XN .NE. 0.0) GO TO 100
C
C        SELECTION OF THE INITIAL APPROXIMATION XN OF X
C                       WHEN A .LT. 1
C
      IF (A .GT. 1.0) GO TO 30
      G = GAMMA(A + 1.0)
      QG = Q*G
      IF (QG .EQ. 0.0) GO TO 560
      IF (QG .GT. 0.6*A) GO TO 20
      B = QG/A
      IF (A .GE. 0.30 .OR. B .LT. 0.35) GO TO 10
         T = EXP(-(B + C))
         U = T*EXP(T)
         XN = T*EXP(U)
         GO TO 100
C
   10 IF (B .GE. 0.45) GO TO 20
      IF (B .EQ. 0.0) GO TO 560
      Y = -ALOG(B)
      S = 0.5 + (0.5 - A)
      T = Y - S*ALOG(Y)
      IF (B .LT. 0.15) GO TO 11
         XN = Y - S*ALOG(T) - ALOG(1.0 + S/(T + 1.0))
         GO TO 100
   11 U = ((T + 2.0*(3.0 - A))*T + (2.0 - A)*(3.0 - A))/
     *          ((T + (5.0 - A))*T + 2.0)
      XN = Y - S*ALOG(T) - ALOG(U)
      GO TO 100
C
   20 U = (P*G)**(1.0/A)
      T = 0.5 + (0.5 - U/(A + 1.0))
      IF (T .GE. E) GO TO 21
         U = EXP(-(C + Q/A))
         XN = U/(1.0 - U/(A + 1.0))
         GO TO 100
   21 XN = U/T
      GO TO 100
C
C        SELECTION OF THE INITIAL APPROXIMATION XN OF X
C                       WHEN A .GT. 1
C
   30 T = Q
      IF (Q .GT. 0.5) T = P
      T = SQRT(-2.0*ALOG(T))
      S = T - (((A3*T + A2)*T + A1)*T + A0)/((((B4*T + B3)*T
     *                + B2)*T + B1)*T + 1.0)
      IF (Q .GT. 0.5) S = -S
C
      RTA = SQRT(A)
      S2 = S*S
      XN = A + S*RTA + (S2 - 1.0)/3.0 + S*(S2 - 7.0)/(36.0*RTA)
     1       - ((3.0*S2 + 7.0)*S2 - 16.0)/(810.0*A)
     2       + S*((9.0*S2 + 256.0)*S2 - 433.0)/(38880.0*A*RTA)
      IF (A .LT. 200.0) GO TO 40
      X = XN
      D = 0.5 + (0.5 - X/A)
      IF (ABS(D) .LE. 2.E-5) RETURN
C
   40 IF (P .LE. 0.5) GO TO 50
      IF (A .GE. 11.0 .OR. XN .LE. 4.0*A) GO TO 200
      B = Q*GAMMA(A)
      IF (B .EQ. 0.0) GO TO 560
      T = A - 1.0
      XN = -ALOG(B) + T*ALOG(XN) - ALOG(1.0 - T/(XN + 1.0))
      GO TO 200
C
   50 IF (XN .GT. 0.11*(A + 1.0)) GO TO 101
      IF (A .GT. 4.0) GO TO 51
      L = INT(A + 0.5)
      IF (P .GT. P0(L)) GO TO 101
   51 W = ALOG(P) + GAMLN(A + 1.0)
      U = EXP(W/A)
      T = SQRT(1.0 - 4.0*U/(A + 1.0))
      X = 2.0*U/(1.0 + T)
      XN = AMAX1(X,XN)
      GO TO 101
C
C                  NEWTON ITERATION USING P
C
  100 IF (P .GT. 0.5) GO TO 200
  101 IF (P .LE. 1.E7*XMIN) GO TO 550
  102 IF (A .LT. AMAX) GO TO 110
      D = 0.5 + (0.5 - XN/A)
      IF (ABS(D) .LE. E2) GO TO 550
C
  110 IF (IERR .GE. 20) GO TO 530
      IERR = IERR + 1
      CALL GRATIO(A,XN,PN,QN,0)
      IF (PN .EQ. 0.0 .OR. QN .EQ. 0.0) GO TO 550
      R = RCOMP(A,XN)
      IF (R .EQ. 0.0) GO TO 550
      X = XN*(1.0 + (P - PN)/R)
      IF (X .LE. 0.0) GO TO 540
C
      XNM1 = XN
      XN = X
      IF (ABS(XN - XNM1) .GT. TOL*X) GO TO 102
      IF (ABS(P - PN) .LE. TOL*P) RETURN
      IF (XN .EQ. XNM1) GO TO 550
      GO TO 102
C
C                  NEWTON ITERATION USING Q
C
  200 IF (Q .LE. 1.E7*XMIN) GO TO 550
  201 IF (A .LT. AMAX) GO TO 210
      D = 0.5 + (0.5 - XN/A)
      IF (ABS(D) .LE. E2) GO TO 550
C
  210 IF (IERR .GE. 20) GO TO 530
      IERR = IERR + 1
      CALL GRATIO(A,XN,PN,QN,0)
      IF (PN .EQ. 0.0 .OR. QN .EQ. 0.0) GO TO 550
      R = RCOMP(A,XN)
      IF (R .EQ. 0.0) GO TO 550
      X = XN*(1.0 + (QN - Q)/R)
      IF (X .LE. 0.0) GO TO 540
C
      XNM1 = XN
      XN = X
      IF (ABS(XN - XNM1) .GT. TOL*X) GO TO 201
      IF (ABS(Q - QN) .LE. TOL*Q) RETURN
      IF (XN .EQ. XNM1) GO TO 550
      GO TO 201
C
C                       SPECIAL CASES
C
  400 X = XMAX
      RETURN
C
  410 X = -ALOG(Q)
      RETURN
C
C                       ERROR RETURN
C
  500 IERR = -2
      RETURN
C
  510 IERR = -4
      RETURN
C
  520 IERR = -5
      RETURN
C
  530 IERR = -6
      RETURN
C
  540 IERR = -7
      RETURN
C
  550 X = XN
      IERR = -8
      RETURN
C
  560 X = XMAX
      IERR = -8
      RETURN
      END
