[gammaresult,x]=gamma(x);

function [gammaresult,x]=gamma(x);
gammaresult=[];
persistent dxrel first firstCall gamma gcs i n ngcs pi sinpiy sq2pil xmax xmin y ; if isempty(firstCall),firstCall=1;end; 

if isempty(dxrel), dxrel=0; end;
if isempty(gammaresult), gammaresult=0; end;
if isempty(gcs), gcs=zeros(1,23); 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(i), i=0; end;
if isempty(n), n=0; end;
if isempty(ngcs), ngcs=0; end;
%***BEGIN PROLOGUE  GAMMA
%***PURPOSE  Compute the complete Gamma function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7A
%***TYPE      SINGLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C)
%***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% GAMMA computes the gamma function at X, where X is not 0, -1, -2, ....
% GAMMA and X are single precision.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  CSEVL, GAMLIM, INITS, R1MACH, 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  GAMMA
if isempty(first), first=false; end;
if firstCall,   gcs(1)=[.008571195590989331e0];  end;
if firstCall,   gcs(2)=[.004415381324841007e0];  end;
if firstCall,   gcs(3)=[.05685043681599363e0];  end;
if firstCall,   gcs(4)=[-.004219835396418561e0];  end;
if firstCall,   gcs(5)=[.001326808181212460e0];  end;
if firstCall,   gcs(6)=[-.0001893024529798880e0];  end;
if firstCall,   gcs(7)=[.0000360692532744124e0];  end;
if firstCall,   gcs(8)=[-.0000060567619044608e0];  end;
if firstCall,   gcs(9)=[.0000010558295463022e0];  end;
if firstCall,   gcs(10)=[-.0000001811967365542e0];  end;
if firstCall,   gcs(11)=[.0000000311772496471e0];  end;
if firstCall,   gcs(12)=[-.0000000053542196390e0];  end;
if firstCall,   gcs(13)=[.0000000009193275519e0];  end;
if firstCall,   gcs(14)=[-.0000000001577941280e0];  end;
if firstCall,   gcs(15)=[.0000000000270798062e0];  end;
if firstCall,   gcs(16)=[-.0000000000046468186e0];  end;
if firstCall,   gcs(17)=[.0000000000007973350e0];  end;
if firstCall,   gcs(18)=[-.0000000000001368078e0];  end;
if firstCall,   gcs(19)=[.0000000000000234731e0];  end;
if firstCall,   gcs(20)=[-.0000000000000040274e0];  end;
if firstCall,   gcs(21)=[.0000000000000006910e0];  end;
if firstCall,   gcs(22)=[-.0000000000000001185e0];  end;
if firstCall,   gcs(23)=[.0000000000000000203e0];  end;
if firstCall,   pi=[3.14159265358979324e0];  end;
% SQ2PIL IS LOG (SQRT (2.*PI) )
if firstCall,   sq2pil=[0.91893853320467274e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%
% LANL DEPENDENT CODE REMOVED 81.02.04
%
%***FIRST EXECUTABLE STATEMENT  GAMMA
if( first )
%
% ---------------------------------------------------------------------
% INITIALIZE.  FIND LEGAL BOUNDS FOR X, AND DETERMINE THE NUMBER OF
% TERMS IN THE SERIES REQUIRED TO ATTAIN AN ACCURACY TEN TIMES BETTER
% THAN MACHINE PRECISION.
%
[ngcs ,gcs]=inits(gcs,23,0.1.*r1mach(3));
%
[xmin,xmax]=gamlim(xmin,xmax);
dxrel = sqrt(r1mach(4));
%
% ---------------------------------------------------------------------
% FINISH INITIALIZATION.  START EVALUATING GAMMA(X).
%
end;
first = false;
%
y = abs(x);
if( y>10.0 )
%
% COMPUTE GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X).
%
if( x>xmax )
xermsg('SLATEC','GAMMA','X SO BIG GAMMA OVERFLOWS',3,2);
end;
%
gammaresult = 0.;
if( x<xmin )
xermsg('SLATEC','GAMMA','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;
%
gammaresult = exp((y-0.5).*log(y)-y+sq2pil+r9lgmc(y));
if( x>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( abs((x-fix(x-0.5))./x)<dxrel )
xermsg('SLATEC','GAMMA','ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER',1,1);
end;
%
sinpiy = sin(pi.*y);
if( sinpiy==0. )
xermsg('SLATEC','GAMMA','X IS A NEGATIVE INTEGER',4,2);
end;
%
gammaresult = -pi./(y.*sinpiy.*gammaresult);
else;
%
% COMPUTE GAMMA(X) FOR ABS(X) .LE. 10.0.  REDUCE INTERVAL AND
% FIND GAMMA(1+Y) FOR 0. .LE. Y .LT. 1. FIRST OF ALL.
%
n = fix(x);
if( x<0. )
n = fix(n - 1);
end;
y = x - n;
n = fix(n - 1);
gammaresult = 0.9375 + csevl(2..*y-1.,gcs,ngcs);
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.
%
for i = 1 : n;
gammaresult =(y+i).*gammaresult;
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.
%
n = fix(-n);
if( x==0. )
xermsg('SLATEC','GAMMA','X IS 0',4,2);
end;
if( x<0. && x+n-2==0. )
xermsg('SLATEC','GAMMA','X IS A NEGATIVE INTEGER',4,2);
end;
if( x<(-0.5) && abs((x-fix(x-0.5))./x)<dxrel )
xermsg('SLATEC','GAMMA','ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',1,1);
end;
%
for i = 1 : n;
gammaresult = gammaresult./(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 GAMR