Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

cgamma_new.m
```% GAMMA  Gamma function valid in the entire complex plane.
%        Accuracy is 15 significant digits along the real axis
%        and 13 significant digits elsewhere.
%        This routine uses a superb Lanczos series
%        approximation for the complex Gamma function.
%
%        z may be complex and of any size.
%        Also  n! = prod(1:n) = gamma(n+1)
%
%usage: [f] = gamma(z)
%
%tested on versions 6.0 and 5.3.1 under Sun Solaris 5.5.1
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931-944
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%            W. J. Cody "An Overview of Software Development for Special
%            Functions", 1975
%
%
%Paul Godfrey
%pgodfrey@intersil.com
%http://winnie.fit.edu/~gabdo/gamma.txt
%Sept 11, 2001

siz = size(z);
z=z(:);
zz=z;

f = 0.*z; % reserve space in advance

p=find(real(z)&lt;0);
if ~isempty(p)
z(p)=-z(p);
end

% 15 sig. digits for 0&lt;=real(z)&lt;=171
% coeffs should sum to about g*g/2+23/24

g=607/128; % best results when 4&lt;=g&lt;=5

c = [  0.99999999999999709182;
57.156235665862923517;
-59.597960355475491248;
14.136097974741747174;
-0.49191381609762019978;
.33994649984811888699e-4;
.46523628927048575665e-4;
-.98374475304879564677e-4;
.15808870322491248884e-3;
-.21026444172410488319e-3;
.21743961811521264320e-3;
-.16431810653676389022e-3;
.84418223983852743293e-4;
-.26190838401581408670e-4;
.36899182659531622704e-5];

%Num Recipes used g=5 with 7 terms
%for a less effective approximation

z=z-1;
zh =z+0.5;
zgh=zh+g;
%trick for avoiding FP overflow above z=141
zp=zgh.^(zh*0.5);

ss=0.0;
for pp=size(c,1)-1:-1:1
ss=ss+c(pp+1)./(z+pp);
end

%sqrt(2Pi)
sq2pi=  2.5066282746310005024157652848110;
f=(sq2pi*(c(1)+ss)).*((zp.*exp(-zgh)).*zp);

f(z==0 | z==1) = 1.0;

if ~isempty(p)
f(p)=-pi./(zz(p).*f(p).*sin(pi*zz(p)));
end

p=find(round(zz)==zz &amp; imag(zz)==0 &amp; real(zz)&lt;=0);
if ~isempty(p)
f(p)=Inf;
end

f=reshape(f,siz);

return

```