dgamit.f

Go to the documentation of this file.
00001 *DECK DGAMIT
00002       DOUBLE PRECISION FUNCTION DGAMIT (A, X)
00003 C***BEGIN PROLOGUE  DGAMIT
00004 C***PURPOSE  Calculate Tricomi's form of the incomplete Gamma function.
00005 C***LIBRARY   SLATEC (FNLIB)
00006 C***CATEGORY  C7E
00007 C***TYPE      DOUBLE PRECISION (GAMIT-S, DGAMIT-D)
00008 C***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
00009 C             SPECIAL FUNCTIONS, TRICOMI
00010 C***AUTHOR  Fullerton, W., (LANL)
00011 C***DESCRIPTION
00012 C
00013 C   Evaluate Tricomi's incomplete Gamma function defined by
00014 C
00015 C   DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
00016 C              T**(A-1.)
00017 C
00018 C   for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
00019 C   GAMMA(X) is the complete gamma function of X.
00020 C
00021 C   DGAMIT is evaluated for arbitrary real values of A and for non-
00022 C   negative values of X (even though DGAMIT is defined for X .LT.
00023 C   0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite,
00024 C   which is a fatal error.
00025 C
00026 C   The function and both arguments are DOUBLE PRECISION.
00027 C
00028 C   A slight deterioration of 2 or 3 digits accuracy will occur when
00029 C   DGAMIT is very large or very small in absolute value, because log-
00030 C   arithmic variables are used.  Also, if the parameter  A  is very
00031 C   close to a negative integer (but not a negative integer), there is
00032 C   a loss of accuracy, which is reported if the result is less than
00033 C   half machine precision.
00034 C
00035 C***REFERENCES  W. Gautschi, A computational procedure for incomplete
00036 C                 gamma functions, ACM Transactions on Mathematical
00037 C                 Software 5, 4 (December 1979), pp. 466-481.
00038 C               W. Gautschi, Incomplete gamma functions, Algorithm 542,
00039 C                 ACM Transactions on Mathematical Software 5, 4
00040 C                 (December 1979), pp. 482-489.
00041 C***ROUTINES CALLED  D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS,
00042 C                    DLNGAM, XERCLR, XERMSG
00043 C***REVISION HISTORY  (YYMMDD)
00044 C   770701  DATE WRITTEN
00045 C   890531  Changed all specific intrinsics to generic.  (WRB)
00046 C   890531  REVISION DATE from Version 3.2
00047 C   891214  Prologue converted to Version 4.0 format.  (BAB)
00048 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
00049 C   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
00050 C***END PROLOGUE  DGAMIT
00051       DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
00052      1  BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, D9LGIT,
00053      2  DLNGAM, D9LGIC
00054       LOGICAL FIRST
00055       SAVE ALNEPS, SQEPS, BOT, FIRST
00056       DATA FIRST /.TRUE./
00057 C***FIRST EXECUTABLE STATEMENT  DGAMIT
00058       IF (FIRST) THEN
00059          ALNEPS = -LOG (D1MACH(3))
00060          SQEPS = SQRT(D1MACH(4))
00061          BOT = LOG (D1MACH(1))
00062       ENDIF
00063       FIRST = .FALSE.
00064 C
00065       IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMIT', 'X IS NEGATIVE'
00066      +   , 2, 2)
00067 C
00068       IF (X.NE.0.D0) ALX = LOG (X)
00069       SGA = 1.0D0
00070       IF (A.NE.0.D0) SGA = SIGN (1.0D0, A)
00071       AINTA = AINT (A + 0.5D0*SGA)
00072       AEPS = A - AINTA
00073 C
00074       IF (X.GT.0.D0) GO TO 20
00075       DGAMIT = 0.0D0
00076       IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
00077       RETURN
00078 C
00079  20   IF (X.GT.1.D0) GO TO 30
00080       IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1,
00081      1  SGNGAM)
00082       DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
00083       RETURN
00084 C
00085  30   IF (A.LT.X) GO TO 40
00086       T = D9LGIT (A, X, DLNGAM(A+1.0D0))
00087       IF (T.LT.BOT) CALL XERCLR
00088       DGAMIT = EXP (T)
00089       RETURN
00090 C
00091  40   ALNG = D9LGIC (A, X, ALX)
00092 C
00093 C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
00094 C
00095       H = 1.0D0
00096       IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50
00097 C
00098       CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM)
00099       T = LOG (ABS(A)) + ALNG - ALGAP1
00100       IF (T.GT.ALNEPS) GO TO 60
00101 C
00102       IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T)
00103       IF (ABS(H).GT.SQEPS) GO TO 50
00104 C
00105       CALL XERCLR
00106       CALL XERMSG ('SLATEC', 'DGAMIT', 'RESULT LT HALF PRECISION', 1,
00107      +   1)
00108 C
00109  50   T = -A*ALX + LOG(ABS(H))
00110       IF (T.LT.BOT) CALL XERCLR
00111       DGAMIT = SIGN (EXP(T), H)
00112       RETURN
00113 C
00114  60   T = T - A*ALX
00115       IF (T.LT.BOT) CALL XERCLR
00116       DGAMIT = -SGA * SGNGAM * EXP(T)
00117       RETURN
00118 C
00119       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines