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

Go to the documentation of this file.
00001       subroutine xgammainc (a, x, result)
00002 
00003 c -- jwe, based on DGAMIT.
00004 c
00005 c -- Do a better job than dgami for large values of x.
00006 
00007       double precision a, x, result
00008       intrinsic exp, log, sqrt, sign, aint
00009       external dgami, dlngam, d9lgit, d9lgic, d9gmit
00010 
00011 C     external dgamr
00012 C     DOUBLE PRECISION DGAMR
00013 
00014       DOUBLE PRECISION AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
00015      $     BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, D9GMIT,
00016      $     D9LGIC, D9LGIT, DLNGAM, DGAMI 
00017 
00018       LOGICAL FIRST
00019 
00020       SAVE ALNEPS, SQEPS, BOT, FIRST
00021 
00022       DATA FIRST /.TRUE./
00023 
00024       if (x .eq. 0.0d0) then
00025 
00026         if (a .eq. 0.0d0) then
00027           result = 1.0d0
00028         else
00029           result = 0.0d0
00030         endif
00031 
00032       else
00033 
00034       IF (FIRST) THEN
00035          ALNEPS = -LOG (D1MACH(3))
00036          SQEPS = SQRT(D1MACH(4))
00037          BOT = LOG (D1MACH(1))
00038       ENDIF
00039       FIRST = .FALSE.
00040 C
00041       IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'XGMAINC', 'X IS NEGATIVE'
00042      +   , 2, 2)
00043 C
00044       IF (X.NE.0.D0) ALX = LOG (X)
00045       SGA = 1.0D0
00046       IF (A.NE.0.D0) SGA = SIGN (1.0D0, A)
00047       AINTA = AINT (A + 0.5D0*SGA)
00048       AEPS = A - AINTA
00049 C
00050 C      IF (X.GT.0.D0) GO TO 20
00051 C      DGAMIT = 0.0D0
00052 C      IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
00053 C      RETURN
00054 C
00055  20   IF (X.GT.1.D0) GO TO 30
00056       IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1,
00057      1  SGNGAM)
00058 C      DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
00059       result = exp (a*alx + log (D9GMIT (A, X, ALGAP1, SGNGAM, ALX)))
00060       RETURN
00061 C
00062  30   IF (A.LT.X) GO TO 40
00063       T = D9LGIT (A, X, DLNGAM(A+1.0D0))
00064       IF (T.LT.BOT) CALL XERCLR
00065 C      DGAMIT = EXP (T)
00066       result = EXP (a*alx + T)
00067       RETURN
00068 C
00069  40   ALNG = D9LGIC (A, X, ALX)
00070 C
00071 C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
00072 C
00073       H = 1.0D0
00074       IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50
00075 C
00076       CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM)
00077       T = LOG (ABS(A)) + ALNG - ALGAP1
00078       IF (T.GT.ALNEPS) GO TO 60
00079 C
00080       IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T)
00081       IF (ABS(H).GT.SQEPS) GO TO 50
00082 C
00083       CALL XERCLR
00084       CALL XERMSG ('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1,
00085      +   1)
00086 C
00087 C 50   T = -A*ALX + LOG(ABS(H))
00088 C      IF (T.LT.BOT) CALL XERCLR
00089 C      DGAMIT = SIGN (EXP(T), H)
00090  50   result = H
00091       RETURN
00092 C
00093 C 60   T = T - A*ALX
00094  60   IF (T.LT.BOT) CALL XERCLR
00095       result = -SGA * SGNGAM * EXP(T)
00096       RETURN
00097 
00098       endif
00099       return
00100       end