/************************************************************************
*                                                                       *
* Purpose:                                                              *
*     randX is a library for random number generation.                  *
*     The integer randXseed should be initialized to an arbitrary       *
* integer prior to the first call to the desired function. The calling  *
* program should not alter the value of randXseed between subsequent    *
* calls.                                                                *
*                                                                       *
* Author:                                                               *
*    Frederico R. B. Cruz                                               *
*    Departamento de Estatistica                                        *
*    Universidade Federal de Minas Gerais                               *
*    Brazil                                                             *
*    E-mail: fcruz@est.ufmg.br                                          *
*                                                                       *
* Version:                                                              *
*    1.00                                                               *
*                                                                       *
* Date:                                                                 *
*    March/2002                                                         *
*                                                                       *
************************************************************************/

#ifndef RANDX_C
#define RANDX_C
#include <math.h>
#define RANDX_PI 3.14159265358979323846
#define RANDX_EE 2.71828182845904523536

/* uniform [0,1]                                                       */
float  rUnif(int *randXseed);
float  rUnif2(unsigned long *randXseed1, unsigned long *randXseed2);

/* normal(mu=0,sd=1)                                                   */
float  rNorm(int *randXseed);

/* exponetial(lambda)                                                  */
float  rExpo(int *randXseed, float lambda);

/* gamma(alpha)                                                        */
double rGamma(int *randXseed, double alpha);

/* beta(alpha,beta)                                                    */
double rBeta(int *randXseed, double alpha, double beta);


/************************************************************************
*                                                                       *
*        The rUnif function is a uniform random number generator based  *
*  on theory and suggestions given in Forsythe, G.E., Malcolm, M.A., &  *
*  Moter, C.B. Computer Methods For Mathematical Computations,          *
*  Prentice-Hall, 1977.                                                 *
*        The integer randXseed should be initialized to an arbitrary    *
* integer prior to the first call to rUnif. The calling program should  *
* not alter the value of randXseed between subsequent calls to rUnif.   *
* Values of rUnif will be returned in the interval (0,1).               *
*                                                                       *
************************************************************************/

float rUnif(int *randXseed) {
    /* initialized data */
    static int m2 = 0;
    static int two = 2;
    /* local variables */
    static int m;
    static float s;
    static float halfm;
    static int a, c, mc;
    /* if first entry then ... */
    if (m2 == 0) {
        /* compute machine integer word length */
	m = 1;
	do {
            m2 = m;
            m = two * m2;
	} while( m > m2 );
	halfm = (float) m2;
        /* compute multiplier and increment for linear  */
        /* congruential method                          */
        a = ( (int) ( halfm * atan(1.) / 8.) << 3 ) + 5;
        c = ( (int) ( halfm * ( 0.5 - sqrt(3.) / 6.)) << 1 ) + 1;
        mc = m2 - c + m2;
        /*   s is the scale factor for converting to floating point */
	s = 0.5 / halfm;
    }
    /* compute next random number */
    *randXseed *= a;
    /*  the following statement is for computers which do not */
    /*  allow integer overflow on addition                    */
    if (*randXseed > mc)
       *randXseed = *randXseed - m2 - m2;
    *randXseed += c;
    /* the following statement is for computers where the word */
    /* length for addition is greater than for multiplication  */
    if (*randXseed / 2 > m2)
       *randXseed = *randXseed - m2 - m2;
    /* the following statement is for computers where integer    */
    /* overflow affects the sign bit                             */
    if (*randXseed < 0)
       *randXseed = *randXseed + m2 + m2;
    return ( (float) (*randXseed) * s );
}

/************************************************************************
*                                                                       *
*        George Marsaglia's uniform random number generator             *
*        The integers randXseed1 and randXseed2 should be initialized   *
* to an arbitrary integer prior to the first call to rUnif2. The        *
* calling program should not alter these values between subsequent      *
* calls to rUnif.                                                       *
*                                                                       *
************************************************************************/

float  rUnif2(unsigned long *randXseed1, unsigned long *randXseed2) {
  static unsigned long s1new, s2new;
  s1new = ((*randXseed1=36969*(*randXseed1&65535)+(*randXseed1>>16))<<16);
  s2new = ((*randXseed2=18000*(*randXseed2&65535)+(*randXseed2>>16))&65535);
  return ((s1new+s2new)*2.32830643708E-10);
}

/************************************************************************
#include <stdio.h>
#include "uni.c"

int main() {
  int replic = 1000;
  int randXseed = 123456;
  unsigned long randXseed1=362436069, randXseed2=521288629;
  int i;
  for (i=0; i<replic; i++){
      fprintf(stdout, "%20.18f\t%20.18f\t%20.18f\n",
              rUnif(&randXseed), rUnif2(&randXseed1,&randXseed2), UNI);
  }
  return 0;
}
************************************************************************/


/************************************************************************
*                                                                       *
*        The rNorm function is a normal random number generator.        *
*        The integer randXseed should be initialized to an arbitrary    *
* integer prior to the first call to rNorm. The calling program should  *
* not alter the value of randXseed between subsequent calls to rNorm.   *
* Values of rNorm will be returned following N(0,1).                    *
*                                                                       *
************************************************************************/

float rNorm(int *randXseed) {
  static float y1, y2;
  while(1) {
      y1=-log(rUnif(randXseed));
      y2=-log(rUnif(randXseed));
      if (y2-(y1-1)*(y1-1)/2 >= 0) {
         if (rUnif(randXseed) > 0.5)
            return -y1;
         else
            return y1;
      }
  }
}

/************************************************************************
#include <stdio.h>
int main() {
 int replic = 100000;
 float mu = 0;
 float sigma = 1;
 int seed = 123456;
 float numb;
 float sum, sum2;
 int i;
 sum = 0.0;
 sum2 = 0.0;
 for (i=0; i<replic; i++){
   numb = mu + rNorm(&seed)*sigma;
   sum += numb;
   sum2 += (numb*numb);
 }
 printf("X ~ N(%f,%f) had E(x)=%f and Var(x)=%f over %d replics.\n",
        mu, sigma*sigma, sum/replic, (sum2-(sum*sum)/replic)/replic, replic);
 return 0;
}
************************************************************************/

/************************************************************************
*                                                                       *
*        The rExpo function is an exponential random number generator.  *
*        The integer randXseed should be initialized to an arbitrary    *
* integer prior to the first call to erand. The calling program should  *
* not alter the value of randXseed between subsequent calls to erand.   *
* Values of erand will be returned following E(lambda).                 *
*                                                                       *
************************************************************************/

float rExpo(int *randXseed, float lambda) {
  return (-log(rUnif(randXseed))/lambda);
}

/************************************************************************
#include <stdio.h>
int main() {
 int replic = 10000;
 float lambda = 2;
 int seed = 123456;
 float numb;
 float sum = 0.0;
 float sum2 = 0.0;
 int i;
 int cont;
 for (i=0; i<replic; i++){
   numb  = rExpo(&seed,lambda);
   sum += numb;
   sum2 += (numb*numb);
 }
 printf("X ~ E(%f) had E(x)=%f^-1 and Var(x)=%f^-2 over %d replics.\n",
        lambda, 1/(sum/replic), 1/sqrt((sum2-(sum*sum)/replic)/replic), replic);
 sum = 0.0;
 cont = 0;
 while (sum <= replic){
   numb  = rExpo(&seed,lambda);
   sum += numb;
   cont++;
 }
 printf("There were %d events over %d time units.\n", cont, replic);
 return 0;
}
************************************************************************/

/************************************************************************
*                                                                       *
*        The rGamma function is a gamma random number generator.        *
*        The integer randXseed should be initialized to an arbitrary    *
* integer prior to the first call to erand. The calling program should  *
* not alter the value of randXseed between subsequent calls to erand.   *
* Values of erand will be returned following G(alpha).                  *
*                                                                       *
************************************************************************/

double rGamma(int *randXseed, double alpha) {
   static double r1,r2,aa,x,w,c1,c2,c3,c4,c5;
   if (alpha<=0.)
      return 0.;
   if (alpha == 1.)
      return rExpo(randXseed,1.);
   if (alpha<1) {
      aa=(alpha+RANDX_EE)/RANDX_EE;
      do {
         r1=rUnif(randXseed);
         r2=rUnif(randXseed);
         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(1);
   }
   else {
      c1=alpha-1;
      c2=(alpha-1./(6.*alpha))/c1;
      c3=2./c1;
      c4=c3+2.;
      c5=1./sqrt(alpha);
      do {
         do {
            r1=rUnif(randXseed);
            r2=rUnif(randXseed);
            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 (1);
   }
}

/************************************************************************
#include <stdio.h>
int main() {
 int replic = 10000;
 int seed = 123456;
 int i;
 double alpha = 1.5;
 for (i=0; i<replic; i++){
   fprintf(stdout, "%f\n", rGamma(&seed,alpha));
 }
 return 0;
}
************************************************************************/

/************************************************************************
*                                                                       *
*        The rBeta function is a beta random number generator.          *
*        The integer randXseed should be initialized to an arbitrary    *
* integer prior to the first call to erand. The calling program should  *
* not alter the value of randXseed between subsequent calls to erand.   *
* Values of erand will be returned following B(alpha,beta).             *
*                                                                       *
************************************************************************/
double rBeta(int *randXseed, double alpha, double beta) {
   double r1;
   if (alpha <=0. || beta <= 0.)
      return 0.;
   r1=rGamma(randXseed,alpha);
   return r1/(r1+rGamma(randXseed,beta));
}

/************************************************************************
#include <stdio.h>
int main() {
 int replic = 10000;
 int seed = 123456;
 int i;
 double alpha = 2.0;
 double beta = 3.0;
 for (i=0; i<replic; i++){
   fprintf(stdout, "%f\n", rBeta(&seed,alpha,beta));
 }
 return 0;
}
************************************************************************/

#endif