#include <stdlib.h>
#include <math.h>	
#define PI 3.141592653589793
#define EE 2.718281828459046
/* seed(int s) - sets seed
   runi() - uniform[0,1]
   runir(a,b) - uniform[a,b]
   runii(na,nb) - uniform on integers na to nb
   rnor(mu,sd) 
   rpois(mu) 
   rexp(lambda) - exponential(lamda) mean = 1/lambda
   rcauchy(location,scale)
   rgamma(alpha)
   rbeta(alpha,beta) 
   rstudent(t) 
   rlogistic(location,scale)
   rlnorm(norm,norsd)
   rbin(n,p)
   rweibull(lambda)
   rf(t,u)
   rchisq(t) */
   
double rxxx1;
int ixxx1=0;
/*************************************/
double runi()
   {
   return ((double)rand() + 0.5)/((double)RAND_MAX + 1.0);
   }
 /* drand not working on DEC for some reason
   {
   double drand48();
   static int idnum = -1;
   return drand48();
   }
 */
/*************************************/
void seed(s)
int s;
   {
   int i;
   for(i=0;i<abs(s);i++) runi();
   /*
   void srand48();
   srand48((long)s);
   */
   }
/*************************************/
void setseedS(s)
int s;
   {
   int i;
   for(i=0;i<abs(s);i++) runi();
   /*
   void srand48();
   srand48((long)s);
   */
   }
/*************************************/
double runir(a,b)
double a,b;
{return (b-a)*runi()+a;}
/*************************************/
int runii(na,nb)
int na,nb;
{return floor((double)(runi()*(nb-na+1.)))+na;}
/*************************************/
double rnor(mu,sd)
double mu,sd;
{
   double e,v1,v2,w;
   if(ixxx1==0){
      do{
         v1=2*runi()-1.;
         v2=2*runi()-1.;
         w=v1*v1+v2*v2;      
      }while(w>1.);
      e=sqrt((-2.*log(w))/w);
      rxxx1=v1*e;
      ixxx1=1;
      return v2*e*sd+mu;
   }
   else{
      ixxx1=0;
      return rxxx1*sd+mu;
   }
}
/*************************************/
int rpois(mu)
double mu;
{ 
   double p=1.;
   int n=0;
   if(mu<=0)return 0;
   mu=exp(-mu);
   do{
      p=p*runi();
      n++;
   }while(p>=mu);
   return n-1;
}
/*************************************/
double rexp(lambda)
double lambda;
{return -log(runi())/lambda;}
/*************************************/
double rcauchy(loc,scale)
double loc,scale;
{return tan((runi()-0.5)*PI)*scale+loc;}
/*************************************/
double rncauchy(loc,scale)
double loc,scale;
{
 double nscale;
 nscale = scale*.67;
 return tan((runi()-0.5)*PI)*nscale+loc;
}
/*************************************/
double rgamma(alpha)
double alpha;
{
   double r1,r2,aa,x,w,c1,c2,c3,c4,c5;
   if(alpha<=0.)return 0.;
   if(alpha == 1.)return rexp(1.);
   if(alpha<1){
      aa=(alpha+EE)/EE;
      do{
         r1=runi();
         r2=runi();
         if(r1>1./aa){
            x = -log(aa*(1.-r1)/alpha);
            if(r2<pow(x,(alpha-1.)))return x;
         }
         else{
            x = pow((aa*r1),(1./alpha));
            if(r2<exp(-x))return x;
         }
      }while(r2<2);
   }
   else{
      c1=alpha-1;
      c2=(alpha-1./(6.*alpha))/c1;
      c3=2./c1;
      c4=c3+2.;
      c5=1./sqrt(alpha);
      do{
         do{
            r1=runi();
            r2=runi();
            if(alpha>2.5)r1=r2+c5*(1.-1.86*r1);
         }while(r1<=0 || r1 >= 1);
         w=c2*r2/r1;
         if(c3*r1+w+1/w <= c4)return c1*w;
         if(c3*log(r1)-log(w)+w<1)return c1*w;
      }while(r2<2);
   }
}
/*************************************/
double rbeta(alpha,beta)
double alpha,beta;
{
   double r1;
   if(alpha <=0. || beta <= 0.) return 0.;
   r1=rgamma(alpha);
   return r1/(r1+rgamma(beta));
}
/*************************************/
double rlogistic(loc,scale)
double loc,scale;
{return 1./(1.-runi())-1;}
/*************************************/
double rlnorm(norm,norsd)
double norm,norsd;
{return exp(rnor(norm,norsd));}
/*************************************/
int rbin(n,p)
int n;
double p;
{
   int j=0,k;
   for(k=0;k<n;k++)if(runi()<p)j++;
   return j;
}
/*************************************/
double rweibull(gamma)
double gamma;
{
   if(gamma<=0)return 0.;
   return pow(rexp(1.),(1/gamma));
}
/*************************************/
double rchisq(t)
double t;
{return rgamma(t/2)*2.;}
/*************************************/
double rf(t,u)
double t,u;
{
   if(t<=0 || u <= 0)return 0.;
   return rchisq(t)*u/(t*rchisq(u));
}
/*************************************/
double rstudent(t)
double t;
{
   double r1,r2;
   if(t<=0.)return 0.;
   return rnor(0.,1.)/sqrt(rchisq(t)/t);
}
/*************************************/