function [dgammaresult,x]=dgamma(x);
dgammaresult=[];
persistent dxrel first firstCall gamcs i n ngam pi sinpiy sq2pil xmax xmin y ; if isempty(firstCall),firstCall=1;end;
;
if isempty(i), i=0; end;
if isempty(n), n=0; end;
if isempty(ngam), ngam=0; end;
%***BEGIN PROLOGUE DGAMMA
%***PURPOSE Compute the complete Gamma function.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C7A
%***TYPE doubleprecision (GAMMA-S, DGAMMA-D, CGAMMA-C)
%***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% DGAMMA(X) calculates the doubleprecision complete Gamma function
% for doubleprecision argument X.
%
% Series for GAM on the interval 0. to 1.00000E+00
% with weighted error 5.79E-32
% log weighted error 31.24
% significant figures required 30.00
% decimal places required 32.05
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG
%***REVISION HISTORY (YYMMDD)
% 770601 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890911 Removed unnecessary intrinsics. (WRB)
% 890911 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)
% 920618 Removed space from variable name. (RWC, WRB)
%***end PROLOGUE DGAMMA
if isempty(gamcs), gamcs=zeros(1,42); end;
if isempty(dxrel), dxrel=0; end;
if isempty(pi), pi=0; end;
if isempty(sinpiy), sinpiy=0; end;
if isempty(sq2pil), sq2pil=0; end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
if isempty(y), y=0; end;
if isempty(first), first=false; end;
%
if firstCall, gamcs(1)=[+.8571195590989331421920062399942d-2]; end;
if firstCall, gamcs(2)=[+.4415381324841006757191315771652d-2]; end;
if firstCall, gamcs(3)=[+.5685043681599363378632664588789d-1]; end;
if firstCall, gamcs(4)=[-.4219835396418560501012500186624d-2]; end;
if firstCall, gamcs(5)=[+.1326808181212460220584006796352d-2]; end;
if firstCall, gamcs(6)=[-.1893024529798880432523947023886d-3]; end;
if firstCall, gamcs(7)=[+.3606925327441245256578082217225d-4]; end;
if firstCall, gamcs(8)=[-.6056761904460864218485548290365d-5]; end;
if firstCall, gamcs(9)=[+.1055829546302283344731823509093d-5]; end;
if firstCall, gamcs(10)=[-.1811967365542384048291855891166d-6]; end;
if firstCall, gamcs(11)=[+.3117724964715322277790254593169d-7]; end;
if firstCall, gamcs(12)=[-.5354219639019687140874081024347d-8]; end;
if firstCall, gamcs(13)=[+.9193275519859588946887786825940d-9]; end;
if firstCall, gamcs(14)=[-.1577941280288339761767423273953d-9]; end;
if firstCall, gamcs(15)=[+.2707980622934954543266540433089d-10]; end;
if firstCall, gamcs(16)=[-.4646818653825730144081661058933d-11]; end;
if firstCall, gamcs(17)=[+.7973350192007419656460767175359d-12]; end;
if firstCall, gamcs(18)=[-.1368078209830916025799499172309d-12]; end;
if firstCall, gamcs(19)=[+.2347319486563800657233471771688d-13]; end;
if firstCall, gamcs(20)=[-.4027432614949066932766570534699d-14]; end;
if firstCall, gamcs(21)=[+.6910051747372100912138336975257d-15]; end;
if firstCall, gamcs(22)=[-.1185584500221992907052387126192d-15]; end;
if firstCall, gamcs(23)=[+.2034148542496373955201026051932d-16]; end;
if firstCall, gamcs(24)=[-.3490054341717405849274012949108d-17]; end;
if firstCall, gamcs(25)=[+.5987993856485305567135051066026d-18]; end;
if firstCall, gamcs(26)=[-.1027378057872228074490069778431d-18]; end;
if firstCall, gamcs(27)=[+.1762702816060529824942759660748d-19]; end;
if firstCall, gamcs(28)=[-.3024320653735306260958772112042d-20]; end;
if firstCall, gamcs(29)=[+.5188914660218397839717833550506d-21]; end;
if firstCall, gamcs(30)=[-.8902770842456576692449251601066d-22]; end;
if firstCall, gamcs(31)=[+.1527474068493342602274596891306d-22]; end;
if firstCall, gamcs(32)=[-.2620731256187362900257328332799d-23]; end;
if firstCall, gamcs(33)=[+.4496464047830538670331046570666d-24]; end;
if firstCall, gamcs(34)=[-.7714712731336877911703901525333d-25]; end;
if firstCall, gamcs(35)=[+.1323635453126044036486572714666d-25]; end;
if firstCall, gamcs(36)=[-.2270999412942928816702313813333d-26]; end;
if firstCall, gamcs(37)=[+.3896418998003991449320816639999d-27]; end;
if firstCall, gamcs(38)=[-.6685198115125953327792127999999d-28]; end;
if firstCall, gamcs(39)=[+.1146998663140024384347613866666d-28]; end;
if firstCall, gamcs(40)=[-.1967938586345134677295103999999d-29]; end;
if firstCall, gamcs(41)=[+.3376448816585338090334890666666d-30]; end;
if firstCall, gamcs(42)=[-.5793070335782135784625493333333d-31]; end;
if firstCall, pi=[3.14159265358979323846264338327950d0]; end;
if firstCall, sq2pil=[0.91893853320467274178032973640562d0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DGAMMA
if( first )
[ngam ,gamcs]=initds(gamcs,42,0.1.*real(d1mach(3)));
%
[xmin,xmax]=dgamlm(xmin,xmax);
dxrel = sqrt(d1mach(4));
end;
first = false;
%
y = abs(x);
if( y>10.0d0 )
%
% GAMMA(X) FOR ABS(X) .GT. 10.0. RECALL Y = ABS(X).
%
if( x>xmax )
xermsg('SLATEC','DGAMMA','X SO BIG GAMMA OVERFLOWS',3,2);
end;
%
dgammaresult = 0.0d0;
if( x<xmin )
xermsg('SLATEC','DGAMMA','X SO SMALL GAMMA UNDERFLOWS',2,1);
end;
if( x<xmin )
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;
%
dgammaresult = exp((y-0.5d0).*log(y)-y+sq2pil+d9lgmc(y));
if( x>0.0d0 )
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;
%
if( abs((x-fix(x-0.5D0))./x)<dxrel )
xermsg('SLATEC','DGAMMA','ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER',1,1);
end;
%
sinpiy = sin(pi.*y);
if( sinpiy==0.0D0 )
xermsg('SLATEC','DGAMMA','X IS A NEGATIVE INTEGER',4,2);
end;
%
dgammaresult = -pi./(y.*sinpiy.*dgammaresult);
else;
%
% COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND. REDUCE INTERVAL AND FIND
% GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL.
%
n = fix(x);
if( x<0.0d0 )
n = fix(n - 1);
end;
y = x - n;
n = fix(n - 1);
dgammaresult = 0.9375d0 + dcsevl(2.0d0.*y-1.0d0,gamcs,ngam);
if( n==0 )
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;
%
if( n>0 )
%
% GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0
%
for i = 1 : n;
dgammaresult =(y+i).*dgammaresult;
end; i = fix(n+1);
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
else;
%
% COMPUTE GAMMA(X) FOR X .LT. 1.0
%
n = fix(-n);
if( x==0.0D0 )
xermsg('SLATEC','DGAMMA','X IS 0',4,2);
end;
if( x<0.0 && x+n-2==0.0D0 )
xermsg('SLATEC','DGAMMA','X IS A NEGATIVE INTEGER',4,2);
end;
if( x<(-0.5D0) && abs((x-fix(x-0.5D0))./x)<dxrel )
xermsg('SLATEC','DGAMMA','ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',1,1);
end;
%
for i = 1 : n;
dgammaresult = dgammaresult./(x+i-1);
end; i = fix(n+1);
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;
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 DGAMR
