/** -----------------------------------------------------------------------
* These are C lang routines for generating random variables from six       
* discrete distributions.                                                   
*                                                                           
* There are generators for the following 6 discrete distributions:          
*                                                                           
*      Generator           Range     Mean         Variance                  
*                                                                           
*      Bernoulli(p)        0..1      p            p*(1-p)                   
*      Binomial(n, p)      0..n      n*p          n*p*(1-p)                 
*      Equilikely(a, b)    a..b      (a+b)/2      ((b-a+1)*(b-a+1)-1)/12    
*      Geometric(p)        0...      p/(1-p)      p/((1-p)*(1-p))           
*      Pascal(n, p)        0...      n*p/(1-p)    n*p/((1-p)*(1-p))         
*      Poisson(m)          0...      m            m                         
*                                                                           
* Name            : DGU.C                                                   
* Purpose         : Discrete Generator Routines                             
* Author          : Steve Park                                              
* Language        : Turbo Pascal, 5.0                                       
* Latest Revision : 09-22-90                                                
* Reference       : Lecture Notes on Simulation, by Steve Park              
* Converted to C  : David W. Geyer  09-02-91                                
* ------------------------------------------------------------------------- 
* NOTE - must link with routine rng.o 
*/

/* include files */
#include <math.h>                          /* for log(), exp() */

/* =======================================================================
*                             long Bernoulli()                              
* Returns 1 with probability p or 0 with probability 1 - p.
*   NOTE: use 0 < p < 1                                      
* ========================================================================= 
*/
long Bernoulli(double p)
{ double Random();

  if (Random() < (1 - p))
    return((long) 0);
  else
    return((long) 1);
}  /* Bernoulli */

/** ========================================================================
*                             long Binomial()                              
*  Returns a binomial distributed integer between 0 and n inclusive.
*   NOTE: use n > 0 and 0 < p < 1                                     
* ========================================================================= 
*/
long Binomial(long n, double p)
{ long i, x = 0, Bernoulli(double);

  for(i=0;i < n;i++)
    x += Bernoulli(p);
  return(x);
}  /* Binomial */

/** ======================================================================
*                            long Equilikely()                             
* Returns an equilikely distributed integer between a and b inclusive.
*   NOTE: use a < b                                                     
* ========================================================================= 
*/
long Equilikely(long a, long b)
{ double Random();

  return(a + (long)((b - a + 1) * Random()));
}  /* Equilikely */

/** ======================================================================= 
*                            long Geometric()                              
* Returns a geometric distributed non-negative integer.
*   NOTE: use 0 < p < 1                                  
* ========================================================================= 

long Geometric(double p)
{ double Random();

  return( (long)( (log(1 - Random()) / log(p)) ) );
}  /* Geometric */

/** =======================================================================
*                            long Pascal()                                  
* Returns a Pascal distributed non-negative integer.
*   NOTE: use n > 0 and 0 < p < 1                      
* ========================================================================= 
*/
long Pascal(long n, double p)
{ long i, x = 0, Geometric(double);

  for (i=0;i < n;i++)
    x += Geometric(p);
  return(x);
}  /* Pascal */

/** =======================================================================
*                            long Poisson()                                 
* Returns a Poisson distributed non-negative integer.
*   NOTE: use m > 0                                     
* ========================================================================= 
*/
long Poisson(double m)
{ long   x = 0;
  double t, Random();

  t = exp(m);
  do {
    x++;
    t *= Random();
  } while (t >= 1);
  return(x - 1);
}  /* Poisson */