| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119 |
- *DECK DGAMIT
- DOUBLE PRECISION FUNCTION DGAMIT (A, X)
- C***BEGIN PROLOGUE DGAMIT
- C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function.
- C***LIBRARY SLATEC (FNLIB)
- C***CATEGORY C7E
- C***TYPE DOUBLE PRECISION (GAMIT-S, DGAMIT-D)
- C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
- C SPECIAL FUNCTIONS, TRICOMI
- C***AUTHOR Fullerton, W., (LANL)
- C***DESCRIPTION
- C
- C Evaluate Tricomi's incomplete Gamma function defined by
- C
- C DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
- C T**(A-1.)
- C
- C for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
- C GAMMA(X) is the complete gamma function of X.
- C
- C DGAMIT is evaluated for arbitrary real values of A and for non-
- C negative values of X (even though DGAMIT is defined for X .LT.
- C 0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite,
- C which is a fatal error.
- C
- C The function and both arguments are DOUBLE PRECISION.
- C
- C A slight deterioration of 2 or 3 digits accuracy will occur when
- C DGAMIT is very large or very small in absolute value, because log-
- C arithmic variables are used. Also, if the parameter A is very
- C close to a negative integer (but not a negative integer), there is
- C a loss of accuracy, which is reported if the result is less than
- C half machine precision.
- C
- C***REFERENCES W. Gautschi, A computational procedure for incomplete
- C gamma functions, ACM Transactions on Mathematical
- C Software 5, 4 (December 1979), pp. 466-481.
- C W. Gautschi, Incomplete gamma functions, Algorithm 542,
- C ACM Transactions on Mathematical Software 5, 4
- C (December 1979), pp. 482-489.
- C***ROUTINES CALLED D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS,
- C DLNGAM, XERCLR, XERMSG
- C***REVISION HISTORY (YYMMDD)
- C 770701 DATE WRITTEN
- C 890531 Changed all specific intrinsics to generic. (WRB)
- C 890531 REVISION DATE from Version 3.2
- C 891214 Prologue converted to Version 4.0 format. (BAB)
- C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
- C 920528 DESCRIPTION and REFERENCES sections revised. (WRB)
- C***END PROLOGUE DGAMIT
- DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
- 1 BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, D9LGIT,
- 2 DLNGAM, D9LGIC
- LOGICAL FIRST
- SAVE ALNEPS, SQEPS, BOT, FIRST
- DATA FIRST /.TRUE./
- C***FIRST EXECUTABLE STATEMENT DGAMIT
- IF (FIRST) THEN
- ALNEPS = -LOG (D1MACH(3))
- SQEPS = SQRT(D1MACH(4))
- BOT = LOG (D1MACH(1))
- ENDIF
- FIRST = .FALSE.
- C
- IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMIT', 'X IS NEGATIVE'
- + , 2, 2)
- C
- IF (X.NE.0.D0) ALX = LOG (X)
- SGA = 1.0D0
- IF (A.NE.0.D0) SGA = SIGN (1.0D0, A)
- AINTA = AINT (A + 0.5D0*SGA)
- AEPS = A - AINTA
- C
- IF (X.GT.0.D0) GO TO 20
- DGAMIT = 0.0D0
- IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
- RETURN
- C
- 20 IF (X.GT.1.D0) GO TO 30
- IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1,
- 1 SGNGAM)
- DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
- RETURN
- C
- 30 IF (A.LT.X) GO TO 40
- T = D9LGIT (A, X, DLNGAM(A+1.0D0))
- IF (T.LT.BOT) CALL XERCLR
- DGAMIT = EXP (T)
- RETURN
- C
- 40 ALNG = D9LGIC (A, X, ALX)
- C
- C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
- C
- H = 1.0D0
- IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50
- C
- CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM)
- T = LOG (ABS(A)) + ALNG - ALGAP1
- IF (T.GT.ALNEPS) GO TO 60
- C
- IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T)
- IF (ABS(H).GT.SQEPS) GO TO 50
- C
- CALL XERCLR
- CALL XERMSG ('SLATEC', 'DGAMIT', 'RESULT LT HALF PRECISION', 1,
- + 1)
- C
- 50 T = -A*ALX + LOG(ABS(H))
- IF (T.LT.BOT) CALL XERCLR
- DGAMIT = SIGN (EXP(T), H)
- RETURN
- C
- 60 T = T - A*ALX
- IF (T.LT.BOT) CALL XERCLR
- DGAMIT = -SGA * SGNGAM * EXP(T)
- RETURN
- C
- END
|
