00001 *DECK GAMLIM 00002 SUBROUTINE GAMLIM (XMIN, XMAX) 00003 C***BEGIN PROLOGUE GAMLIM 00004 C***PURPOSE Compute the minimum and maximum bounds for the argument in 00005 C the Gamma function. 00006 C***LIBRARY SLATEC (FNLIB) 00007 C***CATEGORY C7A, R2 00008 C***TYPE SINGLE PRECISION (GAMLIM-S, DGAMLM-D) 00009 C***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS 00010 C***AUTHOR Fullerton, W., (LANL) 00011 C***DESCRIPTION 00012 C 00013 C Calculate the minimum and maximum legal bounds for X in GAMMA(X). 00014 C XMIN and XMAX are not the only bounds, but they are the only non- 00015 C trivial ones to calculate. 00016 C 00017 C Output Arguments -- 00018 C XMIN minimum legal value of X in GAMMA(X). Any smaller value of 00019 C X might result in underflow. 00020 C XMAX maximum legal value of X in GAMMA(X). Any larger value will 00021 C cause overflow. 00022 C 00023 C***REFERENCES (NONE) 00024 C***ROUTINES CALLED R1MACH, XERMSG 00025 C***REVISION HISTORY (YYMMDD) 00026 C 770401 DATE WRITTEN 00027 C 890531 Changed all specific intrinsics to generic. (WRB) 00028 C 890531 REVISION DATE from Version 3.2 00029 C 891214 Prologue converted to Version 4.0 format. (BAB) 00030 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 00031 C***END PROLOGUE GAMLIM 00032 C***FIRST EXECUTABLE STATEMENT GAMLIM 00033 ALNSML = LOG(R1MACH(1)) 00034 XMIN = -ALNSML 00035 DO 10 I=1,10 00036 XOLD = XMIN 00037 XLN = LOG(XMIN) 00038 XMIN = XMIN - XMIN*((XMIN+0.5)*XLN - XMIN - 0.2258 + ALNSML) 00039 1 / (XMIN*XLN + 0.5) 00040 IF (ABS(XMIN-XOLD).LT.0.005) GO TO 20 00041 10 CONTINUE 00042 CALL XERMSG ('SLATEC', 'GAMLIM', 'UNABLE TO FIND XMIN', 1, 2) 00043 C 00044 20 XMIN = -XMIN + 0.01 00045 C 00046 ALNBIG = LOG(R1MACH(2)) 00047 XMAX = ALNBIG 00048 DO 30 I=1,10 00049 XOLD = XMAX 00050 XLN = LOG(XMAX) 00051 XMAX = XMAX - XMAX*((XMAX-0.5)*XLN - XMAX + 0.9189 - ALNBIG) 00052 1 / (XMAX*XLN - 0.5) 00053 IF (ABS(XMAX-XOLD).LT.0.005) GO TO 40 00054 30 CONTINUE 00055 CALL XERMSG ('SLATEC', 'GAMLIM', 'UNABLE TO FIND XMAX', 2, 2) 00056 C 00057 40 XMAX = XMAX - 0.01 00058 XMIN = MAX (XMIN, -XMAX+1.) 00059 C 00060 RETURN 00061 END
1.7.1
