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

Go to the documentation of this file.
00001 *DECK D9LGIT
00002       DOUBLE PRECISION FUNCTION D9LGIT (A, X, ALGAP1)
00003 C***BEGIN PROLOGUE  D9LGIT
00004 C***SUBSIDIARY
00005 C***PURPOSE  Compute the logarithm of Tricomi's incomplete Gamma
00006 C            function with Perron's continued fraction for large X and
00007 C            A .GE. X.
00008 C***LIBRARY   SLATEC (FNLIB)
00009 C***CATEGORY  C7E
00010 C***TYPE      DOUBLE PRECISION (R9LGIT-S, D9LGIT-D)
00011 C***KEYWORDS  FNLIB, INCOMPLETE GAMMA FUNCTION, LOGARITHM,
00012 C             PERRON'S CONTINUED FRACTION, SPECIAL FUNCTIONS, TRICOMI
00013 C***AUTHOR  Fullerton, W., (LANL)
00014 C***DESCRIPTION
00015 C
00016 C Compute the log of Tricomi's incomplete gamma function with Perron's
00017 C continued fraction for large X and for A .GE. X.
00018 C
00019 C***REFERENCES  (NONE)
00020 C***ROUTINES CALLED  D1MACH, XERMSG
00021 C***REVISION HISTORY  (YYMMDD)
00022 C   770701  DATE WRITTEN
00023 C   890531  Changed all specific intrinsics to generic.  (WRB)
00024 C   890531  REVISION DATE from Version 3.2
00025 C   891214  Prologue converted to Version 4.0 format.  (BAB)
00026 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
00027 C   900720  Routine changed from user-callable to subsidiary.  (WRB)
00028 C***END PROLOGUE  D9LGIT
00029       DOUBLE PRECISION A, X, ALGAP1, AX, A1X, EPS, FK, HSTAR, P, R, S,
00030      1  SQEPS, T, D1MACH
00031       LOGICAL FIRST
00032       SAVE EPS, SQEPS, FIRST
00033       DATA FIRST /.TRUE./
00034 C***FIRST EXECUTABLE STATEMENT  D9LGIT
00035       IF (FIRST) THEN
00036          EPS = 0.5D0*D1MACH(3)
00037          SQEPS = SQRT(D1MACH(4))
00038       ENDIF
00039       FIRST = .FALSE.
00040 C
00041       IF (X .LE. 0.D0 .OR. A .LT. X) CALL XERMSG ('SLATEC', 'D9LGIT',
00042      +   'X SHOULD BE GT 0.0 AND LE A', 2, 2)
00043 C
00044       AX = A + X
00045       A1X = AX + 1.0D0
00046       R = 0.D0
00047       P = 1.D0
00048       S = P
00049       DO 20 K=1,200
00050         FK = K
00051         T = (A+FK)*X*(1.D0+R)
00052         R = T/((AX+FK)*(A1X+FK)-T)
00053         P = R*P
00054         S = S + P
00055         IF (ABS(P).LT.EPS*S) GO TO 30
00056  20   CONTINUE
00057       CALL XERMSG ('SLATEC', 'D9LGIT',
00058      +   'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
00059 C
00060  30   HSTAR = 1.0D0 - X*S/A1X
00061       IF (HSTAR .LT. SQEPS) CALL XERMSG ('SLATEC', 'D9LGIT',
00062      +   'RESULT LESS THAN HALF PRECISION', 1, 1)
00063 C
00064       D9LGIT = -X - ALGAP1 - LOG(HSTAR)
00065       RETURN
00066 C
00067       END