C++生成随机数：Gamma(伽马)分布（gamma distribution）

double cls_random::randomGamma(
    double alpha,
    double lambda)   
{   
    double u, v;   
    double x, y;   
    double b, c;   
    double w, z;   
    bool accept = false;   
    double t;   

    if (alpha > 1.0) 
    {   
        /* Best's rejection algorithm XG for gamma random variates (Best, 1978) */   
        b = alpha - 1;   
        c = 3*alpha-0.75;   
        do {   
            u = cls_random::randomUniform();   
            v = cls_random::randomUniform();   

            w = u*(1-u);   
            y = sqrt(c/w)*(u-0.5);   
            x = b+y;   
            if (x>=0) {   
                z = 64*w*w*w*v*v;
                double zz = 1-2*y*y/x;
                if (z-zz<1e-10)
                {
                    accept = true;
                }
                else
                {
                    accept = false;
                }
                //accept=(z=1-2*y*y/x);   
                //accept=(z==(1-2*y*y/x));   
                if (!accept) 
                {   
                    double logz = log(z);
                    double zzz = 2 * (b*log(x/b)-y);
                    if (logz-zzz<1e-10)
                    {
                        accept = true;
                    }
                    else
                    {
                        accept = false;
                    }
                    //accept=(log(z)==2 * (b*log(x/b)-y));   
                }   
            }   
        } while(!accept);   
        return lambda*x;   
    } 
    else if (alpha == 1.0) 
    {   
        return lambda * cls_random::randomExponential(1/lambda);   
    } 
    else if (alpha < 1.0)
    {   
        double dv = 0.0;
        double b=(exp(1.0)+alpha)/exp(1.0);   
        do
        {   
            double u1=cls_random::randomUniform();
            double p=b*cls_random::randomUniform();  
            double y;   
            if(p>1){   
                y=-log((b-p)/alpha);   
                if(u1<pow(y,alpha-1)){   
                    dv=y;   
                    break;   
                }   
            }else{   
                y=pow(p,1/alpha);   
                if(u1<exp(-y)){   
                    dv=y;   
                    break;   
                }   
            }   
        }while(1);  
        return dv*lambda;
    }   
    return -1;   
}   

期望:E=

方差:V=

wiki: http://zh.wikipedia.org/wiki/%E4%BC%BD%E7%8E%9B%E5%88%86%E5%B8%83