[facresult,n]=fac(n);

function [facresult,n]=fac(n);
facresult=[];
persistent fac facn firstCall nmax sq2pil x xmax xmin ; if isempty(firstCall),firstCall=1;end; 

if isempty(facresult), facresult=0; end;
if isempty(facn), facn=zeros(1,26); 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 isempty(nmax), nmax=0; end;
%***BEGIN PROLOGUE  FAC
%***PURPOSE  Compute the factorial function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C1
%***TYPE      SINGLE PRECISION (FAC-S, DFAC-D)
%***KEYWORDS  FACTORIAL, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% FAC(N) evaluates the factorial function of N.  FAC is single
% precision.  N must be an integer between 0 and 25 inclusive.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  GAMLIM, R9LGMC, 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  FAC
if firstCall,   facn(1)=[1.0e0];  end;
if firstCall,   facn(2)=[1.0e0];  end;
if firstCall,   facn(3)=[2.0e0];  end;
if firstCall,   facn(4)=[6.0e0];  end;
if firstCall,   facn(5)=[24.0e0];  end;
if firstCall,   facn(6)=[120.0e0];  end;
if firstCall,   facn(7)=[720.0e0];  end;
if firstCall,   facn(8)=[5040.0e0];  end;
if firstCall,   facn(9)=[40320.0e0];  end;
if firstCall,   facn(10)=[362880.0e0];  end;
if firstCall,   facn(11)=[3628800.0e0];  end;
if firstCall,   facn(12)=[39916800.0e0];  end;
if firstCall,   facn(13)=[479001600.0e0];  end;
if firstCall,   facn(14)=[6227020800.0e0];  end;
if firstCall,   facn(15)=[87178291200.0e0];  end;
if firstCall,   facn(16)=[1307674368000.0e0];  end;
if firstCall,   facn(17)=[20922789888000.0e0];  end;
if firstCall,   facn(18)=[355687428096000.0e0];  end;
if firstCall,   facn(19)=[6402373705728000.0e0];  end;
if firstCall,   facn(20)=[.12164510040883200e18];  end;
if firstCall,   facn(21)=[.24329020081766400e19];  end;
if firstCall,   facn(22)=[.51090942171709440e20];  end;
if firstCall,   facn(23)=[.11240007277776077e22];  end;
if firstCall,   facn(24)=[.25852016738884977e23];  end;
if firstCall,   facn(25)=[.62044840173323944e24];  end;
if firstCall,   facn(26)=[.15511210043330986e26];  end;
if firstCall,   sq2pil=[0.91893853320467274e0];  end;
if firstCall,   nmax=[0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  FAC
if( nmax==0 )
[xmin,xmax]=gamlim(xmin,xmax);
nmax = fix(xmax - 1.);
end;
%
if( n<0 )
xermsg('SLATEC','FAC','FACTORIAL OF NEGATIVE INTEGER UNDEFINED',1,2);
end;
%
if( n<=25 )
facresult = facn(n+1);
end;
if( n<=25 )
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','FAC','N SO BIG FACTORIAL(N) OVERFLOWS',2,2);
end;
%
x = n + 1;
facresult = exp((x-0.5).*log(x)-x+sq2pil+r9lgmc(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 FC