home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / fortran77_210 / acmtoms / f77 / cgamma404 next >
Encoding:
Text File  |  1992-01-17  |  3.6 KB  |  108 lines
  1. ***********************  COPYRIGHT NOTICE  ********************************
  2. *  The material in this library is copyrighted by the ACM, which grants   *
  3. *  general permission to distribute provided the copies are not made for  *
  4. *  direct commercial advantage.  For details of the copyright and         *
  5. *  dissemination agreement, consult the current issue of TOMS.            *
  6. ***************************************************************************
  7. CDate:           Thu, 17 Oct 91  11:19 BST
  8. CFrom:           CHECAJ@UK.AC.HERIOT-WATT.VAXB
  9. CTo:             KMC@UK.AC.RUTHERFORD.DEC-E
  10. CSubject:        computes the Complex gamma function
  11. C     ALGORITHM 404 COLLECTED ALGORITHMS FROM ACM.
  12. C     ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 01,
  13. C     P. 048.
  14.       FUNCTION CGAMMA(Z)
  15.       COMPLEX Z,ZM,T,TT,SUM,TERM,DEN,CGAMMA,PI,A
  16.       DIMENSION C(12)
  17.       LOGICAL REFLEK
  18. C SET IOUT FOR PROPER OUTPUT CHANNEL OF COMPUTER SYSTEM FOR
  19. C ERROR MESSAGES
  20.       IOUT = 3
  21.       PI = (3.141593,0.0)
  22.       X = REAL(Z)
  23.       Y = AIMAG(Z)
  24. C TOL = LIMIT OF PRECISION OF COMPUTER SYSTEM IN SINGLE PRECISI
  25.       TOL = 1.0E-7
  26.       REFLEK = .TRUE.
  27. C DETERMINE WHETHER Z IS TOO CLOSE TO A POLE
  28. C CHECK WHETHER TOO CLOSE TO ORIGIN
  29.       IF(X.GE.TOL) GO TO 20
  30. C FIND THE NEAREST POLE AND COMPUTE DISTANCE TO IT
  31.       XDIST = X-INT(X-.5)
  32.       ZM = CMPLX(XDIST,Y)
  33.       IF(CABS(ZM).GE.TOL) GO TO 10
  34. C IF Z IS TOO CLOSE TO A POLE, PRINT ERROR MESSAGE AND RETURN
  35. C WITH CGAMMA = (1.E7,0.0E0)
  36.       WRITE(IOUT,900) Z
  37.       CGAMMA = (1.E7,0.E0)
  38.       RETURN
  39. C FOR REAL(Z) NEGATIVE EMPLOY THE REFLECTION FORMULA
  40. C GAMMA(Z) = PI/(SIN(PI*Z)*GAMMA(1-Z))
  41. C AND COMPUTE GAMMA(1-Z).  NOTE REFLEK IS A TAG TO INDICATE THA
  42. C THIS RELATION MUST BE USED LATER.
  43. 10    IF(X.GE.0.0) GO TO 20
  44.       REFLEK = .FALSE.
  45.       Z = (1.0,0.0)-Z
  46.       X = 1.0-X
  47.       Y = -Y
  48. C IF Z IS NOT TOO CLOSE TO A POLE, MAKE REAL(Z)>10 AND ARG(Z)<P
  49. 20    M = 0
  50. 40    IF(X.GE.10.) GO TO 50
  51.       X = X + 1.0
  52.       M = M + 1
  53.       GO TO 40
  54. 50    IF(ABS(Y).LT.X) GO TO 60
  55.       X = X + 1.0
  56.       M = M + 1
  57.       GO TO 50
  58. 60    T = CMPLX(X,Y)
  59.       TT = T*T
  60.       DEN = T
  61. C COEFFICIENTS IN STIRLING*S APPROXIMATION FOR LN(GAMMA(T))
  62.       C(1) = 1./12.
  63.       C(2) = -1./360.
  64.       C(3) = 1./1260.
  65.       C(4) = -1./1680.
  66.       C(5) = 1./1188.
  67.       C(6) = -691./360360.
  68.       C(7) = 1./156.
  69.       C(8) = -3617./122400.
  70.       C(9) = 43867./244188.
  71.       C(10) = -174611./125400.
  72.       C(11) = 77683./5796.
  73.       SUM = (T-(.5,0.0))*CLOG(T)-T+CMPLX(.5*ALOG(2.*3.14159),0.0)
  74.       J = 1
  75. 70    TERM = C(J)/DEN
  76. C TEST REAL AND IMAGINARY PARTS OF LN(GAMMA(Z)) SEPARATELY FOR
  77. C CONVERGENCE.  IF Z IS REAL SKIP IMAGINARY PART OF CHECK.
  78.       IF(ABS(REAL(TERM)/REAL(SUM)).GE.TOL) GO TO 80
  79.       IF(Y.EQ.0.0) GO TO 100
  80.       IF(ABS(AIMAG(TERM)/AIMAG(SUM)).LT.TOL) GO TO 100
  81. 80    SUM = SUM + TERM
  82.       J = J + 1
  83.       DEN = DEN*TT
  84. C TEST FOR NONCONVERGENCE
  85.       IF(J.EQ.12) GO TO 90
  86.       GO TO 70
  87. C STIRLING*S SERIES DID NOT CONVERGE.  PRINT ERROR MESSAGE AND
  88. C PROCEDE.
  89. 90    WRITE(IOUT,910) Z
  90. C RECURSION RELATION USED TO OBTAIN LN(GAMMA(Z))
  91. C LN(GAMMA(Z)) = LN(GAMMA(Z+M)/(Z*(Z+1)*...*(Z+M-1)))
  92. C = LN(GAMMA(Z+M)-LN(Z)-LN(Z+1)-...-LN(Z+M
  93. 100   IF(M.EQ.0) GO TO 120
  94.       DO 110 I = 1,M
  95.       A = CMPLX(I*1.-1.,0.0)
  96. 110   SUM = SUM-CLOG(Z+A)
  97. C CHECK TO SEE IF REFLECTION FORMULA SHOULD BE USED
  98. 120   IF(REFLEK) GO TO 130
  99.       SUM = CLOG(PI/CSIN(PI*Z))-SUM
  100.       Z = (1.0,0.0) -Z
  101. 130   CGAMMA = CEXP(SUM)
  102.       RETURN
  103. 900   FORMAT(1X,2E14.7,10X,49HARGUMENT OF GAMMA FUNCTION IS TOO CLOSE TO
  104.      1 A POLE)
  105. 910   FORMAT(44H ERROR - STIRLING*S SERIES HAS NOT CONVERGED/14X,4HZ = ,
  106.      12E14.7)
  107.       END
  108.