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

Go to the documentation of this file.
00001 *DECK R9GMIT
00002       FUNCTION R9GMIT (A, X, ALGAP1, SGNGAM, ALX)
00003 C***BEGIN PROLOGUE  R9GMIT
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      SINGLE 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  ALNGAM, R1MACH, XERMSG
00019 C***REVISION HISTORY  (YYMMDD)
00020 C   770701  DATE WRITTEN
00021 C   890531  Changed all specific intrinsics to generic.  (WRB)
00022 C   890531  REVISION DATE from Version 3.2
00023 C   891214  Prologue converted to Version 4.0 format.  (BAB)
00024 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
00025 C   900720  Routine changed from user-callable to subsidiary.  (WRB)
00026 C***END PROLOGUE  R9GMIT
00027       SAVE EPS, BOT
00028       DATA EPS, BOT / 2*0.0 /
00029 C***FIRST EXECUTABLE STATEMENT  R9GMIT
00030       IF (EPS.EQ.0.0) EPS = 0.5*R1MACH(3)
00031       IF (BOT.EQ.0.0) BOT = LOG(R1MACH(1))
00032 C
00033       IF (X .LE. 0.0) CALL XERMSG ('SLATEC', 'R9GMIT',
00034      +   'X SHOULD BE GT 0', 1, 2)
00035 C
00036       MA = A + 0.5
00037       IF (A.LT.0.0) MA = A - 0.5
00038       AEPS = A - MA
00039 C
00040       AE = A
00041       IF (A.LT.(-0.5)) AE = AEPS
00042 C
00043       T = 1.0
00044       TE = AE
00045       S = T
00046       DO 20 K=1,200
00047         FK = K
00048         TE = -X*TE/FK
00049         T = TE/(AE+FK)
00050         S = S + T
00051         IF (ABS(T).LT.EPS*ABS(S)) GO TO 30
00052  20   CONTINUE
00053       CALL XERMSG ('SLATEC', 'R9GMIT',
00054      +   'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
00055 C
00056  30   IF (A.GE.(-0.5)) ALGS = -ALGAP1 + LOG(S)
00057       IF (A.GE.(-0.5)) GO TO 60
00058 C
00059       ALGS = -ALNGAM(1.0+AEPS) + LOG(S)
00060       S = 1.0
00061       M = -MA - 1
00062       IF (M.EQ.0) GO TO 50
00063       T = 1.0
00064       DO 40 K=1,M
00065         T = X*T/(AEPS-M-1+K)
00066         S = S + T
00067         IF (ABS(T).LT.EPS*ABS(S)) GO TO 50
00068  40   CONTINUE
00069 C
00070  50   R9GMIT = 0.0
00071       ALGS = -MA*LOG(X) + ALGS
00072       IF (S.EQ.0.0 .OR. AEPS.EQ.0.0) GO TO 60
00073 C
00074       SGNG2 = SGNGAM*SIGN(1.0,S)
00075       ALG2 = -X - ALGAP1 + LOG(ABS(S))
00076 C
00077       IF (ALG2.GT.BOT) R9GMIT = SGNG2*EXP(ALG2)
00078       IF (ALGS.GT.BOT) R9GMIT = R9GMIT + EXP(ALGS)
00079       RETURN
00080 C
00081  60   R9GMIT = EXP(ALGS)
00082       RETURN
00083 C
00084       END