src/numerical/slatec/fortran/dgamlm.f

   1 *DECK DGAMLM
   2       SUBROUTINE DGAMLM (XMIN, XMAX)
   3 C***BEGIN PROLOGUE  DGAMLM
   4 C***PURPOSE  Compute the minimum and maximum bounds for the argument in
   5 C            the Gamma function.
   6 C***LIBRARY   SLATEC (FNLIB)
   7 C***CATEGORY  C7A, R2
   8 C***TYPE      DOUBLE PRECISION (GAMLIM-S, DGAMLM-D)
   9 C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS
  10 C***AUTHOR  Fullerton, W., (LANL)
  11 C***DESCRIPTION
  12 C
  13 C Calculate the minimum and maximum legal bounds for X in gamma(X).
  14 C XMIN and XMAX are not the only bounds, but they are the only non-
  15 C trivial ones to calculate.
  16 C
  17 C             Output Arguments --
  18 C XMIN   double precision minimum legal value of X in gamma(X).  Any
  19 C        smaller value of X might result in underflow.
  20 C XMAX   double precision maximum legal value of X in gamma(X).  Any
  21 C        larger value of X might cause overflow.
  22 C
  23 C***REFERENCES  (NONE)
  24 C***ROUTINES CALLED  D1MACH, XERMSG
  25 C***REVISION HISTORY  (YYMMDD)
  26 C   770601  DATE WRITTEN
  27 C   890531  Changed all specific intrinsics to generic.  (WRB)
  28 C   890531  REVISION DATE from Version 3.2
  29 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  30 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  31 C***END PROLOGUE  DGAMLM
  32       DOUBLE PRECISION XMIN, XMAX, ALNBIG, ALNSML, XLN, XOLD, D1MACH
  33 C***FIRST EXECUTABLE STATEMENT  DGAMLM
  34       ALNSML = LOG(D1MACH(1))
  35       XMIN = -ALNSML
  36       DO 10 I=1,10
  37         XOLD = XMIN
  38         XLN = LOG(XMIN)
  39         XMIN = XMIN - XMIN*((XMIN+0.5D0)*XLN - XMIN - 0.2258D0 + ALNSML)
  40      1    / (XMIN*XLN+0.5D0)
  41         IF (ABS(XMIN-XOLD).LT.0.005D0) GO TO 20
  42  10   CONTINUE
  43       CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMIN', 1, 2)
  44 C
  45  20   XMIN = -XMIN + 0.01D0
  46 C
  47       ALNBIG = LOG (D1MACH(2))
  48       XMAX = ALNBIG
  49       DO 30 I=1,10
  50         XOLD = XMAX
  51         XLN = LOG(XMAX)
  52         XMAX = XMAX - XMAX*((XMAX-0.5D0)*XLN - XMAX + 0.9189D0 - ALNBIG)
  53      1    / (XMAX*XLN-0.5D0)
  54         IF (ABS(XMAX-XOLD).LT.0.005D0) GO TO 40
  55  30   CONTINUE
  56       CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMAX', 2, 2)
  57 C
  58  40   XMAX = XMAX - 0.01D0
  59       XMIN = MAX (XMIN, -XMAX+1.D0)
  60 C
  61       RETURN
  62       END