function gamma=gammarv(alpha,beta,n);
% gammarv(alpha, beta,n) generates a vector of n random numbers
% from a Gamma(alpha,beta) distribution:
%   f(x)=exp(-x/beta)*x^(alpha-1)/[Gamma(alpha)*beta^alpha]   if x>=0
%   f(x)=0 if x<0
% where Gamma(z) is the gamma function.
% Both alpha and beta must be positive.
% Mean(x)=alpha*beta; Var(x)=alpha*beta^2
% See G.S. Fishman, 1973, Concepts and Methods in Discrete Event
% Digital Simulation, Wiley, New York, pp. 208-209.

gamma=[];
for i=1:n
    X=0;
    k=floor(alpha);
    g=alpha-k;
    if k>0
        X=-log(prod(rand(1,k)));
    end;
    if g==0
        gamma=[gamma beta*X];
    else
        a=g;
        b=1-g;
        y=1;
        z=1;
        while y+z>1
           y=rand^(1/a);
           z=rand^(1/b);
        end;
        Y=y/(y+z);      % Y is a Beta(g,1-g) r.v.
        Z=-log(rand);   % Z is a Gamma(1,1) r.v.
        gamma=[gamma beta*(X+Y*Z)];
    end;
end;