ガンマ乱数発生のプログラム

#include <oxstd.h>
main(){
decl alpha,i,u1,u2,x1,x2,a,b,c,x,y;
 // set the value of alpha
alpha=0.5;
// ranseed(100);  // change seed 
for(i=0;i<50;++i)
{

if(alpha > 1.0){
 a=((alpha-1)/exp(1))^((alpha-1)/2);
 b=((alpha+1)/exp(1))^((alpha+1)/2);
  do{
     u1=ranu(1,1);u2=ranu(1,1);
     x1=u1*a;x2=u2*b;
     y=x2/x1;
    }while(2*log(x1)-(alpha-1)*log(y)+y >1);
 }
else if(alpha==1.0){
     u1=ranu(1,1);y=-log(u1);
         }
else{
  do{
     u1=ranu(1,1);u2=ranu(1,1);
     c=1+alpha/exp(1);
     if(u1 <=1/c){
       x=(c*u1)^(1/alpha);y=exp(-x);
           }
     else{
       x=-log(c*(1-u1)/alpha);y=x^(alpha-1);
           }
        } while(u2 >=y);
 }
println(y);
}
}