/***************************************************************************
 * Copyright (C) 1995, 1996, Jun Sun, jsun@junsun.net.                     *
 *                                                                         *
 ***************************************************************************
 */
/* implementation of class Random_Variate */
#include "random_variates.h"
#include "LCG.h"
#include <stdio.h>
#include <stdlib.h>
#include "cppstring.h"

/* const class names */
const String 	RANDOM_VARIATE_CLASS("Random_Variate class");
const String	DISCRETE_RANDOM_VARIATE_CLASS("Discrete_Random_Variate class");
const String	DUNIFORM_CLASS("Duniform class");
const String	BERNOULLI_CLASS("Bernoulli class");
const String	BINOMIAL_CLASS("Binomial class");
const String	GEOMETRIC_CLASS("Geometric class");
const String	NEGBIN_CLASS("Negbin class");
const String	POISSON_CLASS("Poisson class");
const String	CONTINUOUS_RANDOM_VARIATE_CLASS("Continuous_Random_Variate class");
const String	CUNIFORM_CLASS("CUniform class");
const String	EXPONENTIAL_CLASS("Exponential class");
const String	GAMMA_CLASS("Gamma class");
const String	WEIBULL_CLASS("Weibull class");
const String	NORMAL_CLASS("Normal class");
const String	LOGNORMAL_CLASS("Lognormal class");


LCG	sys_rand;
// the system random number generator, default generator for all variates
// in the system.
// the reason for all variates sharing the same generator comes from the 
// fact that the default seed is system time in second, which will vary
// very little for different copies, therefore leads to correlation
// between variates in the system.  On the other hand system time is nice
// to be default seed, so that the results can differ runs from runs.

Random_Variate::Random_Variate()
{ randp = &sys_rand;		// default value 
}


Random_Variate::~Random_Variate()
{ }

// currently termination is the only way to deal with exception
// hopefully in the future, mythrow/catch mechanism can be used
void Random_Variate::mythrow(char* mesg)
{ fprintf(stderr,"%s",mesg);
  exit(-1);
}


void Random_Variate::set_reproducible(int seed)
{ sys_rand.set_seed(seed);}



/********************************************************************/
/*                  helper functions defined here                   */
/********************************************************************/
double factorial(int n)
{ double f = 1.0;
  int i;
  for(i=2;i<=n;i++)
    f *= (double)i;
  return f;
}

double combination(int m, int n)  // m choose n
{ if (n>m) return -1.0;
  else
    return factorial(m)/factorial(n)/factorial(m-n);
}