GNU Octave: libcruft/slatec-fn/d9gmit.f Source File

Go to the documentation of this file.
00001 *DECK D9GMIT
00002       DOUBLE PRECISION FUNCTION D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
00003 C***BEGIN PROLOGUE  D9GMIT
00004 C***SUBSIDIARY
00005 C***PURPOSE  Compute Tricomi's incomplete Gamma function for small
00006 C            arguments.
00007 C***LIBRARY   SLATEC (FNLIB)
00008 C***CATEGORY  C7E
00009 C***TYPE      DOUBLE PRECISION (R9GMIT-S, D9GMIT-D)
00010 C***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, SMALL X,
00011 C             SPECIAL FUNCTIONS, TRICOMI
00012 C***AUTHOR  Fullerton, W., (LANL)
00013 C***DESCRIPTION
00014 C
00015 C Compute Tricomi's incomplete gamma function for small X.
00016 C
00017 C***REFERENCES  (NONE)
00018 C***ROUTINES CALLED  D1MACH, DLNGAM, XERMSG
00019 C***REVISION HISTORY  (YYMMDD)
00020 C   770701  DATE WRITTEN
00021 C   890531  Changed all specific intrinsics to generic.  (WRB)
00022 C   890911  Removed unnecessary intrinsics.  (WRB)
00023 C   890911  REVISION DATE from Version 3.2
00024 C   891214  Prologue converted to Version 4.0 format.  (BAB)
00025 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
00026 C   900720  Routine changed from user-callable to subsidiary.  (WRB)
00027 C***END PROLOGUE  D9GMIT
00028       DOUBLE PRECISION A, X, ALGAP1, SGNGAM, ALX, AE, AEPS, ALGS, ALG2,
00029      1  BOT, EPS, FK, S, SGNG2, T, TE, D1MACH, DLNGAM
00030       LOGICAL FIRST
00031       SAVE EPS, BOT, FIRST
00032       DATA FIRST /.TRUE./
00033 C***FIRST EXECUTABLE STATEMENT  D9GMIT
00034       IF (FIRST) THEN
00035          EPS = 0.5D0*D1MACH(3)
00036          BOT = LOG (D1MACH(1))
00037       ENDIF
00038       FIRST = .FALSE.
00039 C
00040       IF (X .LE. 0.D0) CALL XERMSG ('SLATEC', 'D9GMIT',
00041      +   'X SHOULD BE GT 0', 1, 2)
00042 C
00043       MA = A + 0.5D0
00044       IF (A.LT.0.D0) MA = A - 0.5D0
00045       AEPS = A - MA
00046 C
00047       AE = A
00048       IF (A.LT.(-0.5D0)) AE = AEPS
00049 C
00050       T = 1.D0
00051       TE = AE
00052       S = T
00053       DO 20 K=1,200
00054         FK = K
00055         TE = -X*TE/FK
00056         T = TE/(AE+FK)
00057         S = S + T
00058         IF (ABS(T).LT.EPS*ABS(S)) GO TO 30
00059  20   CONTINUE
00060       CALL XERMSG ('SLATEC', 'D9GMIT',
00061      +   'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
00062 C
00063  30   IF (A.GE.(-0.5D0)) ALGS = -ALGAP1 + LOG(S)
00064       IF (A.GE.(-0.5D0)) GO TO 60
00065 C
00066       ALGS = -DLNGAM(1.D0+AEPS) + LOG(S)
00067       S = 1.0D0
00068       M = -MA - 1
00069       IF (M.EQ.0) GO TO 50
00070       T = 1.0D0
00071       DO 40 K=1,M
00072         T = X*T/(AEPS-(M+1-K))
00073         S = S + T
00074         IF (ABS(T).LT.EPS*ABS(S)) GO TO 50
00075  40   CONTINUE
00076 C
00077  50   D9GMIT = 0.0D0
00078       ALGS = -MA*LOG(X) + ALGS
00079       IF (S.EQ.0.D0 .OR. AEPS.EQ.0.D0) GO TO 60
00080 C
00081       SGNG2 = SGNGAM * SIGN (1.0D0, S)
00082       ALG2 = -X - ALGAP1 + LOG(ABS(S))
00083 C
00084       IF (ALG2.GT.BOT) D9GMIT = SGNG2 * EXP(ALG2)
00085       IF (ALGS.GT.BOT) D9GMIT = D9GMIT + EXP(ALGS)
00086       RETURN
00087 C
00088  60   D9GMIT = EXP (ALGS)
00089       RETURN
00090 C
00091       END