! REV PROCESSED 213 LINES OF CODE. PROGRAM DONE.
REAL FUNCTION GAMMA1(X)
use da_control, only : pi
!D DOUBLE PRECISION FUNCTION DGAMMA(X)
!----------------------------------------------------------------------
!
! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
! MACHINE-DEPENDENT CONSTANTS.
!
!
!*******************************************************************
!*******************************************************************
!
! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
!
! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
! GAMMA(XBIG) = BETA**MAXEXP
! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
! APPROXIMATELY BETA**MAXEXP
! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
! 1.0+EPS .GT. 1.0
! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
! 1/XMININ IS MACHINE REPRESENTABLE
!
! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
!
! BETA MAXEXP XBIG
!
! CRAY-1 (S.P.) 2 8191 966.961
! CYBER 180/855
! UNDER NOS (S.P.) 2 1070 177.803
! IEEE (IBM/XT,
! SUN, ETC.) (S.P.) 2 128 35.040
! IEEE (IBM/XT,
! SUN, ETC.) (D.P.) 2 1024 171.624
! IBM 3033 (D.P.) 16 63 57.574
! VAX D-FORMAT (D.P.) 2 127 34.844
! VAX G-FORMAT (D.P.) 2 1023 171.489
!
! XINF EPS XMININ
!
! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
! CYBER 180/855
! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
! IEEE (IBM/XT,
! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
! IEEE (IBM/XT,
! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
!
!*******************************************************************
!*******************************************************************
!
! ERROR RETURNS
!
! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
! TO BE FREE OF UNDERFLOW AND OVERFLOW.
!
!
! INTRINSIC FUNCTIONS REQUIRED ARE:
!
! INT, DBLE, EXP, LOG, REAL, SIN
!
!
! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
! (ED.), SPRINGER VERLAG, BERLIN, 1976.
!
! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
! SONS, NEW YORK, 1968.
!
! LATEST MODIFICATION: OCTOBER 12, 1989
!
! AUTHORS: W. J. CODY AND L. STOLTZ
! APPLIED MATHEMATICS DIVISION
! ARGONNE NATIONAL LABORATORY
! ARGONNE, IL 60439
!
!----------------------------------------------------------------------
INTEGER I,N
LOGICAL PARITY
!D DOUBLE PRECISION
REAL :: &
C,CONV,EPS,FACT,HALF,ONE,P,Q,RES,SQRTPI,SUM,TWELVE, &
TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
DIMENSION C(7),P(8),Q(8)
!----------------------------------------------------------------------
! MATHEMATICAL CONSTANTS
!----------------------------------------------------------------------
DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/, &
SQRTPI/0.9189385332046727417803297E0/
!D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
!D 1 SQRTPI/0.9189385332046727417803297D0/,
!----------------------------------------------------------------------
! MACHINE DEPENDENT PARAMETERS
!----------------------------------------------------------------------
DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/, &
XINF/3.4E38/
!D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
!D 1 XINF/1.79D308/
!----------------------------------------------------------------------
! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
! APPROXIMATION OVER (1,2).
!----------------------------------------------------------------------
DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
-3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
-3.61444134186911729807069E+4,6.64561438202405440627855E+4/
DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
-1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
-1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
!D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
!D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
!D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
!D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
!D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
!D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
!D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3,
!D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
!----------------------------------------------------------------------
! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
!----------------------------------------------------------------------
DATA C/-1.910444077728E-03,8.4171387781295E-04, &
-5.952379913043012E-04,7.93650793500350248E-04, &
-2.777777777777681622553E-03,8.333333333333333331554247E-02, &
5.7083835261E-03/
!D DATA C/-1.910444077728D-03,8.4171387781295D-04,
!D 1 -5.952379913043012D-04,7.93650793500350248D-04,
!D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02,
!D 3 5.7083835261D-03/
!----------------------------------------------------------------------
! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
!----------------------------------------------------------------------
CONV(I) = REAL(I)
!D CONV(I) = DBLE(I)
PARITY=.FALSE.
FACT=ONE
N=0
Y=X
IF(Y.LE.ZERO)THEN
!----------------------------------------------------------------------
! ARGUMENT IS NEGATIVE
!----------------------------------------------------------------------
Y=-X
Y1=AINT(Y)
RES=Y-Y1
IF(RES.NE.ZERO)THEN
IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
FACT=-PI/SIN(PI*RES)
Y=Y+ONE
ELSE
RES=XINF
GOTO 900
ENDIF
ENDIF
!----------------------------------------------------------------------
! ARGUMENT IS POSITIVE
!----------------------------------------------------------------------
IF(Y.LT.EPS)THEN
!----------------------------------------------------------------------
! ARGUMENT .LT. EPS
!----------------------------------------------------------------------
IF(Y.GE.XMININ)THEN
RES=ONE/Y
ELSE
RES=XINF
GOTO 900
ENDIF
ELSEIF(Y.LT.TWELVE)THEN
Y1=Y
IF(Y.LT.ONE)THEN
!----------------------------------------------------------------------
! 0.0 .LT. ARGUMENT .LT. 1.0
!----------------------------------------------------------------------
Z=Y
Y=Y+ONE
ELSE
!----------------------------------------------------------------------
! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
!----------------------------------------------------------------------
N=INT(Y)-1
Y=Y-CONV(N)
Z=Y-ONE
ENDIF
!----------------------------------------------------------------------
! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
!----------------------------------------------------------------------
XNUM=ZERO
XDEN=ONE
DO 260 I=1,8
XNUM=(XNUM+P(I))*Z
XDEN=XDEN*Z+Q(I)
260 CONTINUE
RES=XNUM/XDEN+ONE
IF(Y1.LT.Y)THEN
!----------------------------------------------------------------------
! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
!----------------------------------------------------------------------
RES=RES/Y1
ELSEIF(Y1.GT.Y)THEN
!----------------------------------------------------------------------
! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
!----------------------------------------------------------------------
DO 290 I=1,N
RES=RES*Y
Y=Y+ONE
290 CONTINUE
ENDIF
ELSE
!----------------------------------------------------------------------
! EVALUATE FOR ARGUMENT .GE. 12.0,
!----------------------------------------------------------------------
IF(Y.LE.XBIG)THEN
YSQ=Y*Y
SUM=C(7)
DO 350 I=1,6
SUM=SUM/YSQ+C(I)
350 CONTINUE
SUM=SUM/Y-Y+SQRTPI
SUM=SUM+(Y-HALF)*LOG(Y)
RES=EXP(SUM)
ELSE
RES=XINF
GOTO 900
ENDIF
ENDIF
!----------------------------------------------------------------------
! FINAL ADJUSTMENTS AND RETURN
!----------------------------------------------------------------------
IF(PARITY)RES=-RES
IF(FACT.NE.ONE)RES=FACT/RES
900 GAMMA1=RES
!D900 DGAMMA = RES
! ---------- LAST LINE OF GAMMA ----------
END FUNCTION GAMMA1