C     ALGORITHM 404 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 01,
C     P. 048.
      FUNCTION CGAMMA(Z)
      COMPLEX Z,ZM,T,TT,SUM,TERM,DEN,CGAMMA,PI,A
      DIMENSION C(12)
      LOGICAL REFLEK
C SET IOUT FOR PROPER OUTPUT CHANNEL OF COMPUTER SYSTEM FOR
C ERROR MESSAGES
      IOUT = 3
      PI = (3.141593,0.0)
      X = REAL(Z)
      Y = AIMAG(Z)
C TOL = LIMIT OF PRECISION OF COMPUTER SYSTEM IN SINGLE PRECISI
      TOL = 1.0E-7
      REFLEK = .TRUE.
C DETERMINE WHETHER Z IS TOO CLOSE TO A POLE
C CHECK WHETHER TOO CLOSE TO ORIGIN
      IF(X.GE.TOL) GO TO 20
C FIND THE NEAREST POLE AND COMPUTE DISTANCE TO IT
      XDIST = X-INT(X-.5)
      ZM = CMPLX(XDIST,Y)
      IF(CABS(ZM).GE.TOL) GO TO 10
C IF Z IS TOO CLOSE TO A POLE, PRINT ERROR MESSAGE AND RETURN
C WITH CGAMMA = (1.E7,0.0E0)
      WRITE(IOUT,900) Z
      CGAMMA = (1.E7,0.E0)
      RETURN
C FOR REAL(Z) NEGATIVE EMPLOY THE REFLECTION FORMULA
C GAMMA(Z) = PI/(SIN(PI*Z)*GAMMA(1-Z))
C AND COMPUTE GAMMA(1-Z).  NOTE REFLEK IS A TAG TO INDICATE THA
C THIS RELATION MUST BE USED LATER.
10    IF(X.GE.0.0) GO TO 20
      REFLEK = .FALSE.
      Z = (1.0,0.0)-Z
      X = 1.0-X
      Y = -Y
C IF Z IS NOT TOO CLOSE TO A POLE, MAKE REAL(Z)>10 AND ARG(Z)<P
20    M = 0
40    IF(X.GE.10.) GO TO 50
      X = X + 1.0
      M = M + 1
      GO TO 40
50    IF(ABS(Y).LT.X) GO TO 60
      X = X + 1.0
      M = M + 1
      GO TO 50
60    T = CMPLX(X,Y)
      TT = T*T
      DEN = T
C COEFFICIENTS IN STIRLING*S APPROXIMATION FOR LN(GAMMA(T))
      C(1) = 1./12.
      C(2) = -1./360.
      C(3) = 1./1260.
      C(4) = -1./1680.
      C(5) = 1./1188.
      C(6) = -691./360360.
      C(7) = 1./156.
      C(8) = -3617./122400.
      C(9) = 43867./244188.
      C(10) = -174611./125400.
      C(11) = 77683./5796.
      SUM = (T-(.5,0.0))*CLOG(T)-T+CMPLX(.5*ALOG(2.*3.14159),0.0)
      J = 1
70    TERM = C(J)/DEN
C TEST REAL AND IMAGINARY PARTS OF LN(GAMMA(Z)) SEPARATELY FOR
C CONVERGENCE.  IF Z IS REAL SKIP IMAGINARY PART OF CHECK.
      IF(ABS(REAL(TERM)/REAL(SUM)).GE.TOL) GO TO 80
      IF(Y.EQ.0.0) GO TO 100
      IF(ABS(AIMAG(TERM)/AIMAG(SUM)).LT.TOL) GO TO 100
80    SUM = SUM + TERM
      J = J + 1
      DEN = DEN*TT
C TEST FOR NONCONVERGENCE
      IF(J.EQ.12) GO TO 90
      GO TO 70
C STIRLING*S SERIES DID NOT CONVERGE.  PRINT ERROR MESSAGE AND
C PROCEDE.
90    WRITE(IOUT,910) Z
C RECURSION RELATION USED TO OBTAIN LN(GAMMA(Z))
C LN(GAMMA(Z)) = LN(GAMMA(Z+M)/(Z*(Z+1)*...*(Z+M-1)))
C = LN(GAMMA(Z+M)-LN(Z)-LN(Z+1)-...-LN(Z+M
100   IF(M.EQ.0) GO TO 120
      DO 110 I = 1,M
      A = CMPLX(I*1.-1.,0.0)
110   SUM = SUM-CLOG(Z+A)
C CHECK TO SEE IF REFLECTION FORMULA SHOULD BE USED
120   IF(REFLEK) GO TO 130
      SUM = CLOG(PI/CSIN(PI*Z))-SUM
      Z = (1.0,0.0) -Z
130   CGAMMA = CEXP(SUM)
      RETURN
900   FORMAT(1X,2E14.7,10X,49HARGUMENT OF GAMMA FUNCTION IS TOO CLOSE TO
     1 A POLE)
910   FORMAT(44H ERROR - STIRLING*S SERIES HAS NOT CONVERGED/14X,4HZ = ,
     12E14.7)
      END