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

Go to the documentation of this file.
00001 *DECK DGAMLM
00002       SUBROUTINE DGAMLM (XMIN, XMAX)
00003 C***BEGIN PROLOGUE  DGAMLM
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      DOUBLE 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   double precision minimum legal value of X in gamma(X).  Any
00019 C        smaller value of X might result in underflow.
00020 C XMAX   double precision maximum legal value of X in gamma(X).  Any
00021 C        larger value of X might cause overflow.
00022 C
00023 C***REFERENCES  (NONE)
00024 C***ROUTINES CALLED  D1MACH, XERMSG
00025 C***REVISION HISTORY  (YYMMDD)
00026 C   770601  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  DGAMLM
00032       DOUBLE PRECISION XMIN, XMAX, ALNBIG, ALNSML, XLN, XOLD, D1MACH
00033 C***FIRST EXECUTABLE STATEMENT  DGAMLM
00034       ALNSML = LOG(D1MACH(1))
00035       XMIN = -ALNSML
00036       DO 10 I=1,10
00037         XOLD = XMIN
00038         XLN = LOG(XMIN)
00039         XMIN = XMIN - XMIN*((XMIN+0.5D0)*XLN - XMIN - 0.2258D0 + ALNSML)
00040      1    / (XMIN*XLN+0.5D0)
00041         IF (ABS(XMIN-XOLD).LT.0.005D0) GO TO 20
00042  10   CONTINUE
00043       CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMIN', 1, 2)
00044 C
00045  20   XMIN = -XMIN + 0.01D0
00046 C
00047       ALNBIG = LOG (D1MACH(2))
00048       XMAX = ALNBIG
00049       DO 30 I=1,10
00050         XOLD = XMAX
00051         XLN = LOG(XMAX)
00052         XMAX = XMAX - XMAX*((XMAX-0.5D0)*XLN - XMAX + 0.9189D0 - ALNBIG)
00053      1    / (XMAX*XLN-0.5D0)
00054         IF (ABS(XMAX-XOLD).LT.0.005D0) GO TO 40
00055  30   CONTINUE
00056       CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMAX', 2, 2)
00057 C
00058  40   XMAX = XMAX - 0.01D0
00059       XMIN = MAX (XMIN, -XMAX+1.D0)
00060 C
00061       RETURN
00062       END