[dgamrresult,x]=dgamr(x);

function [dgamrresult,x]=dgamr(x);
dgamrresult=[];
persistent alngx irold sgngx ; 

;
if isempty(irold), irold=0; end;
%***BEGIN PROLOGUE  DGAMR
%***PURPOSE  Compute the reciprocal of the Gamma function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7A
%***TYPE      doubleprecision (GAMR-S, DGAMR-D, CGAMR-C)
%***KEYWORDS  FNLIB, RECIPROCAL GAMMA FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% DGAMR(X) calculates the doubleprecision reciprocal of the
% complete Gamma function for doubleprecision argument X.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  DGAMMA, DLGAMS, XERCLR, XGETF, XSETF
%***REVISION HISTORY  (YYMMDD)
%   770701  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)
%   900727  Added EXTERNAL statement.  (WRB)
%***end PROLOGUE  DGAMR
if isempty(alngx), alngx=0; end;
if isempty(sgngx), sgngx=0; end;
%***FIRST EXECUTABLE STATEMENT  DGAMR
dgamrresult = 0.0d0;
if( x<=0.0d0 && fix(x)==x )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
[irold]=xgetf(irold);
xsetf(1);
if( abs(x)>10.0d0 )
%
[x,alngx,sgngx]=dlgams(x,alngx,sgngx);
xerclr;
[irold]=xsetf(irold);
dgamrresult = sgngx.*exp(-alngx);
else;
dgamrresult = 1.0d0./dgamma(x);
xerclr;
[irold]=xsetf(irold);
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DGAMRN