var/da/da_tools/gamma1.f90

   1 ! REV PROCESSED 213 LINES OF CODE. PROGRAM DONE.
   2 REAL FUNCTION GAMMA1(X)
   3    use da_control, only : pi
   4 !D    DOUBLE PRECISION FUNCTION DGAMMA(X)
   5 !----------------------------------------------------------------------
   6 !
   7 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
   8 !   COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
   9 !   THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
  10 !   FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS.  COEFFICIENTS
  11 !   FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
  12 !   THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
  13 !   THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
  14 !   COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
  15 !   MACHINE-DEPENDENT CONSTANTS.
  16 !
  17 !
  18 !*******************************************************************
  19 !*******************************************************************
  20 !
  21 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
  22 !
  23 ! BETA   - RADIX FOR THE FLOATING-POINT REPRESENTATION
  24 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
  25 ! XBIG   - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
  26 !          IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
  27 !                  GAMMA(XBIG) = BETA**MAXEXP
  28 ! XINF   - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
  29 !          APPROXIMATELY BETA**MAXEXP
  30 ! EPS    - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
  31 !          1.0+EPS .GT. 1.0
  32 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
  33 !          1/XMININ IS MACHINE REPRESENTABLE
  34 !
  35 !     APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
  36 !
  37 !                            BETA       MAXEXP        XBIG
  38 !
  39 ! CRAY-1         (S.P.)        2         8191        966.961
  40 ! CYBER 180/855
  41 !   UNDER NOS    (S.P.)        2         1070        177.803
  42 ! IEEE (IBM/XT,
  43 !   SUN, ETC.)   (S.P.)        2          128        35.040
  44 ! IEEE (IBM/XT,
  45 !   SUN, ETC.)   (D.P.)        2         1024        171.624
  46 ! IBM 3033       (D.P.)       16           63        57.574
  47 ! VAX D-FORMAT   (D.P.)        2          127        34.844
  48 ! VAX G-FORMAT   (D.P.)        2         1023        171.489
  49 !
  50 !                            XINF         EPS        XMININ
  51 !
  52 ! CRAY-1         (S.P.)   5.45E+2465   7.11E-15    1.84E-2466
  53 ! CYBER 180/855
  54 !   UNDER NOS    (S.P.)   1.26E+322    3.55E-15    3.14E-294
  55 ! IEEE (IBM/XT,
  56 !   SUN, ETC.)   (S.P.)   3.40E+38     1.19E-7     1.18E-38
  57 ! IEEE (IBM/XT,
  58 !   SUN, ETC.)   (D.P.)   1.79D+308    2.22D-16    2.23D-308
  59 ! IBM 3033       (D.P.)   7.23D+75     2.22D-16    1.39D-76
  60 ! VAX D-FORMAT   (D.P.)   1.70D+38     1.39D-17    5.88D-39
  61 ! VAX G-FORMAT   (D.P.)   8.98D+307    1.11D-16    1.12D-308
  62 !
  63 !*******************************************************************
  64 !*******************************************************************
  65 !
  66 ! ERROR RETURNS
  67 !
  68 !  THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
  69 !     WHEN OVERFLOW WOULD OCCUR.  THE COMPUTATION IS BELIEVED
  70 !     TO BE FREE OF UNDERFLOW AND OVERFLOW.
  71 !
  72 !
  73 !  INTRINSIC FUNCTIONS REQUIRED ARE:
  74 !
  75 !     INT, DBLE, EXP, LOG, REAL, SIN
  76 !
  77 !
  78 ! REFERENCES:  AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
  79 !              FUNCTIONS   W. J. CODY, LECTURE NOTES IN MATHEMATICS,
  80 !              506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
  81 !              (ED.), SPRINGER VERLAG, BERLIN, 1976.
  82 !
  83 !              COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
  84 !              SONS, NEW YORK, 1968.
  85 !
  86 !  LATEST MODIFICATION: OCTOBER 12, 1989
  87 !
  88 !  AUTHORS: W. J. CODY AND L. STOLTZ
  89 !           APPLIED MATHEMATICS DIVISION
  90 !           ARGONNE NATIONAL LABORATORY
  91 !           ARGONNE, IL 60439
  92 !
  93 !----------------------------------------------------------------------
  94       INTEGER I,N
  95       LOGICAL PARITY
  96 !D    DOUBLE PRECISION
  97       REAL :: &
  98          C,CONV,EPS,FACT,HALF,ONE,P,Q,RES,SQRTPI,SUM,TWELVE,  &
  99          TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
 100       DIMENSION C(7),P(8),Q(8)
 101 !----------------------------------------------------------------------
 102 !  MATHEMATICAL CONSTANTS
 103 !----------------------------------------------------------------------
 104       DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/,  &
 105            SQRTPI/0.9189385332046727417803297E0/
 106 !D    DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
 107 !D   1     SQRTPI/0.9189385332046727417803297D0/,
 108 !----------------------------------------------------------------------
 109 !  MACHINE DEPENDENT PARAMETERS
 110 !----------------------------------------------------------------------
 111       DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,  &
 112           XINF/3.4E38/
 113 !D    DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
 114 !D   1     XINF/1.79D308/
 115 !----------------------------------------------------------------------
 116 !  NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
 117 !     APPROXIMATION OVER (1,2).
 118 !----------------------------------------------------------------------
 119       DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
 120              -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
 121              8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
 122              -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
 123       DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
 124             -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
 125               2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
 126             -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
 127 !D    DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
 128 !D   1       -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
 129 !D   2       8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
 130 !D   3       -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
 131 !D    DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
 132 !D   1      -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
 133 !D   2        2.25381184209801510330112D+4,4.75584627752788110767815D+3,
 134 !D   3      -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
 135 !----------------------------------------------------------------------
 136 !  COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
 137 !----------------------------------------------------------------------
 138       DATA C/-1.910444077728E-03,8.4171387781295E-04, &
 139            -5.952379913043012E-04,7.93650793500350248E-04, &
 140            -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
 141             5.7083835261E-03/
 142 !D    DATA C/-1.910444077728D-03,8.4171387781295D-04,
 143 !D   1     -5.952379913043012D-04,7.93650793500350248D-04,
 144 !D   2     -2.777777777777681622553D-03,8.333333333333333331554247D-02,
 145 !D   3      5.7083835261D-03/
 146 !----------------------------------------------------------------------
 147 !  STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
 148 !----------------------------------------------------------------------
 149       CONV(I) = REAL(I)
 150 !D    CONV(I) = DBLE(I)
 151       PARITY=.FALSE.
 152       FACT=ONE
 153       N=0
 154       Y=X
 155       IF(Y.LE.ZERO)THEN
 156 !----------------------------------------------------------------------
 157 !  ARGUMENT IS NEGATIVE
 158 !----------------------------------------------------------------------
 159         Y=-X
 160         Y1=AINT(Y)
 161         RES=Y-Y1
 162         IF(RES.NE.ZERO)THEN
 163           IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
 164           FACT=-PI/SIN(PI*RES)
 165           Y=Y+ONE
 166         ELSE
 167           RES=XINF
 168           GOTO 900
 169         ENDIF
 170       ENDIF
 171 !----------------------------------------------------------------------
 172 !  ARGUMENT IS POSITIVE
 173 !----------------------------------------------------------------------
 174       IF(Y.LT.EPS)THEN
 175 !----------------------------------------------------------------------
 176 !  ARGUMENT .LT. EPS
 177 !----------------------------------------------------------------------
 178         IF(Y.GE.XMININ)THEN
 179           RES=ONE/Y
 180         ELSE
 181           RES=XINF
 182           GOTO 900
 183         ENDIF
 184       ELSEIF(Y.LT.TWELVE)THEN
 185         Y1=Y
 186         IF(Y.LT.ONE)THEN
 187 !----------------------------------------------------------------------
 188 !  0.0 .LT. ARGUMENT .LT. 1.0
 189 !----------------------------------------------------------------------
 190           Z=Y
 191           Y=Y+ONE
 192         ELSE
 193 !----------------------------------------------------------------------
 194 !  1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
 195 !----------------------------------------------------------------------
 196           N=INT(Y)-1
 197           Y=Y-CONV(N)
 198           Z=Y-ONE
 199         ENDIF
 200 !----------------------------------------------------------------------
 201 !  EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
 202 !----------------------------------------------------------------------
 203         XNUM=ZERO
 204         XDEN=ONE
 205         DO 260 I=1,8
 206           XNUM=(XNUM+P(I))*Z
 207           XDEN=XDEN*Z+Q(I)
 208   260   CONTINUE
 209         RES=XNUM/XDEN+ONE
 210         IF(Y1.LT.Y)THEN
 211 !----------------------------------------------------------------------
 212 !  ADJUST RESULT FOR CASE  0.0 .LT. ARGUMENT .LT. 1.0
 213 !----------------------------------------------------------------------
 214           RES=RES/Y1
 215         ELSEIF(Y1.GT.Y)THEN
 216 !----------------------------------------------------------------------
 217 !  ADJUST RESULT FOR CASE  2.0 .LT. ARGUMENT .LT. 12.0
 218 !----------------------------------------------------------------------
 219           DO 290 I=1,N
 220             RES=RES*Y
 221             Y=Y+ONE
 222   290     CONTINUE
 223         ENDIF
 224       ELSE
 225 !----------------------------------------------------------------------
 226 !  EVALUATE FOR ARGUMENT .GE. 12.0,
 227 !----------------------------------------------------------------------
 228         IF(Y.LE.XBIG)THEN
 229           YSQ=Y*Y
 230           SUM=C(7)
 231           DO 350 I=1,6
 232             SUM=SUM/YSQ+C(I)
 233   350     CONTINUE
 234           SUM=SUM/Y-Y+SQRTPI
 235           SUM=SUM+(Y-HALF)*LOG(Y)
 236           RES=EXP(SUM)
 237         ELSE
 238           RES=XINF
 239           GOTO 900
 240         ENDIF
 241       ENDIF
 242 !----------------------------------------------------------------------
 243 !  FINAL ADJUSTMENTS AND RETURN
 244 !----------------------------------------------------------------------
 245       IF(PARITY)RES=-RES
 246       IF(FACT.NE.ONE)RES=FACT/RES
 247   900 GAMMA1=RES
 248 !D900 DGAMMA = RES
 249 ! ---------- LAST LINE OF GAMMA ----------
 250 END FUNCTION GAMMA1
 251