function [dfacresult,n]=dfac(n);
dfacresult=[];
persistent facn firstCall nmax sq2pil x xmax xmin ; if isempty(firstCall),firstCall=1;end;
;
if isempty(nmax), nmax=0; end;
%***BEGIN PROLOGUE DFAC
%***PURPOSE Compute the factorial function.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C1
%***TYPE doubleprecision (FAC-S, DFAC-D)
%***KEYWORDS FACTORIAL, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% DFAC(N) calculates the doubleprecision factorial for integer
% argument N.
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D9LGMC, DGAMLM, XERMSG
%***REVISION HISTORY (YYMMDD)
% 770601 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 DFAC
if isempty(facn), facn=zeros(1,31); end;
if isempty(sq2pil), sq2pil=0; end;
if isempty(x), x=0; end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
if firstCall, facn(1)=[+.100000000000000000000000000000000d+1]; end;
if firstCall, facn(2)=[+.100000000000000000000000000000000d+1]; end;
if firstCall, facn(3)=[+.200000000000000000000000000000000d+1]; end;
if firstCall, facn(4)=[+.600000000000000000000000000000000d+1]; end;
if firstCall, facn(5)=[+.240000000000000000000000000000000d+2]; end;
if firstCall, facn(6)=[+.120000000000000000000000000000000d+3]; end;
if firstCall, facn(7)=[+.720000000000000000000000000000000d+3]; end;
if firstCall, facn(8)=[+.504000000000000000000000000000000d+4]; end;
if firstCall, facn(9)=[+.403200000000000000000000000000000d+5]; end;
if firstCall, facn(10)=[+.362880000000000000000000000000000d+6]; end;
if firstCall, facn(11)=[+.362880000000000000000000000000000d+7]; end;
if firstCall, facn(12)=[+.399168000000000000000000000000000d+8]; end;
if firstCall, facn(13)=[+.479001600000000000000000000000000d+9]; end;
if firstCall, facn(14)=[+.622702080000000000000000000000000d+10]; end;
if firstCall, facn(15)=[+.871782912000000000000000000000000d+11]; end;
if firstCall, facn(16)=[+.130767436800000000000000000000000d+13]; end;
if firstCall, facn(17)=[+.209227898880000000000000000000000d+14]; end;
if firstCall, facn(18)=[+.355687428096000000000000000000000d+15]; end;
if firstCall, facn(19)=[+.640237370572800000000000000000000d+16]; end;
if firstCall, facn(20)=[+.121645100408832000000000000000000d+18]; end;
if firstCall, facn(21)=[+.243290200817664000000000000000000d+19]; end;
if firstCall, facn(22)=[+.510909421717094400000000000000000d+20]; end;
if firstCall, facn(23)=[+.112400072777760768000000000000000d+22]; end;
if firstCall, facn(24)=[+.258520167388849766400000000000000d+23]; end;
if firstCall, facn(25)=[+.620448401733239439360000000000000d+24]; end;
if firstCall, facn(26)=[+.155112100433309859840000000000000d+26]; end;
if firstCall, facn(27)=[+.403291461126605635584000000000000d+27]; end;
if firstCall, facn(28)=[+.108888694504183521607680000000000d+29]; end;
if firstCall, facn(29)=[+.304888344611713860501504000000000d+30]; end;
if firstCall, facn(30)=[+.884176199373970195454361600000000d+31]; end;
if firstCall, facn(31)=[+.265252859812191058636308480000000d+33]; end;
if firstCall, sq2pil=[0.91893853320467274178032973640562d0]; end;
if firstCall, nmax=[0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DFAC
if( nmax==0 )
[xmin,xmax]=dgamlm(xmin,xmax);
nmax = fix(xmax - 1.0d0);
end;
%
if( n<0 )
xermsg('SLATEC','DFAC','FACTORIAL OF NEGATIVE INTEGER UNDEFINED',1,2);
end;
%
if( n<=30 )
dfacresult = facn(n+1);
end;
if( n<=30 )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
if( n>nmax )
xermsg('SLATEC','DFAC','N SO BIG FACTORIAL(N) OVERFLOWS',2,2);
end;
%
x = n + 1;
dfacresult = exp((x-0.5d0).*log(x)-x+sq2pil+d9lgmc(x));
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DFC
