00001 *DECK GAMIT 00002 REAL FUNCTION GAMIT (A, X) 00003 C***BEGIN PROLOGUE GAMIT 00004 C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function. 00005 C***LIBRARY SLATEC (FNLIB) 00006 C***CATEGORY C7E 00007 C***TYPE SINGLE 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 GAMIT = 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 GAMIT is evaluated for arbitrary real values of A and for non- 00022 C negative values of X (even though GAMIT is defined for X .LT. 00023 C 0.0), except that for X = 0 and A .LE. 0.0, GAMIT is infinite, 00024 C which is a fatal error. 00025 C 00026 C The function and both arguments are REAL. 00027 C 00028 C A slight deterioration of 2 or 3 digits accuracy will occur when 00029 C GAMIT 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 ALGAMS, ALNGAM, GAMR, R1MACH, R9GMIT, R9LGIC, 00042 C R9LGIT, 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 GAMIT 00051 LOGICAL FIRST 00052 SAVE ALNEPS, SQEPS, BOT, FIRST 00053 DATA FIRST /.TRUE./ 00054 C***FIRST EXECUTABLE STATEMENT GAMIT 00055 IF (FIRST) THEN 00056 ALNEPS = -LOG(R1MACH(3)) 00057 SQEPS = SQRT(R1MACH(4)) 00058 BOT = LOG(R1MACH(1)) 00059 ENDIF 00060 FIRST = .FALSE. 00061 C 00062 IF (X .LT. 0.0) CALL XERMSG ('SLATEC', 'GAMIT', 'X IS NEGATIVE', 00063 + 2, 2) 00064 C 00065 IF (X.NE.0.0) ALX = LOG(X) 00066 SGA = 1.0 00067 IF (A.NE.0.0) SGA = SIGN (1.0, A) 00068 AINTA = AINT (A+0.5*SGA) 00069 AEPS = A - AINTA 00070 C 00071 IF (X.GT.0.0) GO TO 20 00072 GAMIT = 0.0 00073 IF (AINTA.GT.0.0 .OR. AEPS.NE.0.0) GAMIT = GAMR(A+1.0) 00074 RETURN 00075 C 00076 20 IF (X.GT.1.0) GO TO 40 00077 IF (A.GE.(-0.5) .OR. AEPS.NE.0.0) CALL ALGAMS (A+1.0, ALGAP1, 00078 1 SGNGAM) 00079 GAMIT = R9GMIT (A, X, ALGAP1, SGNGAM, ALX) 00080 RETURN 00081 C 00082 40 IF (A.LT.X) GO TO 50 00083 T = R9LGIT (A, X, ALNGAM(A+1.0)) 00084 IF (T.LT.BOT) CALL XERCLR 00085 GAMIT = EXP(T) 00086 RETURN 00087 C 00088 50 ALNG = R9LGIC (A, X, ALX) 00089 C 00090 C EVALUATE GAMIT IN TERMS OF LOG(GAMIC(A,X)) 00091 C 00092 H = 1.0 00093 IF (AEPS.EQ.0.0 .AND. AINTA.LE.0.0) GO TO 60 00094 CALL ALGAMS (A+1.0, ALGAP1, SGNGAM) 00095 T = LOG(ABS(A)) + ALNG - ALGAP1 00096 IF (T.GT.ALNEPS) GO TO 70 00097 IF (T.GT.(-ALNEPS)) H = 1.0 - SGA*SGNGAM*EXP(T) 00098 IF (ABS(H).GT.SQEPS) GO TO 60 00099 CALL XERCLR 00100 CALL XERMSG ('SLATEC', 'GAMIT', 'RESULT LT HALF PRECISION', 1, 1) 00101 C 00102 60 T = -A*ALX + LOG(ABS(H)) 00103 IF (T.LT.BOT) CALL XERCLR 00104 GAMIT = SIGN (EXP(T), H) 00105 RETURN 00106 C 00107 70 T = T - A*ALX 00108 IF (T.LT.BOT) CALL XERCLR 00109 GAMIT = -SGA*SGNGAM*EXP(T) 00110 RETURN 00111 C 00112 END
1.7.1
