[xmin,xmax]=gamlim(xmin,xmax);

function [xmin,xmax]=gamlim(xmin,xmax);
persistent alnbig alnsml gt i xln xold ; 

if isempty(alnbig), alnbig=0; end;
if isempty(alnsml), alnsml=0; end;
if isempty(xln), xln=0; end;
if isempty(xold), xold=0; end;
if isempty(i), i=0; end;
if isempty(gt), gt=0; end;
%***BEGIN PROLOGUE  GAMLIM
%***PURPOSE  Compute the minimum and maximum bounds for the argument in
%            the Gamma function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7A, R2
%***TYPE      SINGLE PRECISION (GAMLIM-S, DGAMLM-D)
%***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Calculate the minimum and maximum legal bounds for X in GAMMA(X).
% XMIN and XMAX are not the only bounds, but they are the only non-
% trivial ones to calculate.
%
%             Output Arguments --
% XMIN   minimum legal value of X in GAMMA(X).  Any smaller value of
%        X might result in underflow.
% XMAX   maximum legal value of X in GAMMA(X).  Any larger value will
%        cause overflow.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770401  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%***end PROLOGUE  GAMLIM
%***FIRST EXECUTABLE STATEMENT  GAMLIM
alnsml = log(r1mach(1));
xmin = -alnsml;
gt=0;
for i = 1 : 10;
xold = xmin;
xln = log(xmin);
xmin = xmin - xmin.*((xmin+0.5).*xln-xmin-0.2258+alnsml)./(xmin.*xln+0.5);
if( abs(xmin-xold)<0.005 )
gt=1;
break;
end;
end;
if(gt==0)
xermsg('SLATEC','GAMLIM','UNABLE TO FIND XMIN',1,2);
end;
%
xmin = -xmin + 0.01;
%
alnbig = log(r1mach(2));
xmax = alnbig;
gt=0;
for i = 1 : 10;
xold = xmax;
xln = log(xmax);
xmax = xmax - xmax.*((xmax-0.5).*xln-xmax+0.9189-alnbig)./(xmax.*xln-0.5);
if( abs(xmax-xold)<0.005 )
gt=1;
break;
end;
end;
if(gt==0)
xermsg('SLATEC','GAMLIM','UNABLE TO FIND XMAX',2,2);
end;
%
xmax = xmax - 0.01;
xmin = max(xmin,-xmax+1.);
%
end %subroutine gamlim
%DECK GAMLN