#include <stdio.h>
#include <stdlib.h>
//#include <conio.h>
#include <math.h>
/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */

#if (defined (__BORLANDC__) || defined (_MSC_VER)) && ! defined(_WINDOWS_)                   
  #define _GETCH_DEFINED_
#endif


void FatalError(char * ErrorText);


const double SHAT1 = 2.943035529371538573;    // 8/e
const double SHAT2 = 0.8989161620588987408;   // 3-sqrt(12/e)

static const int FAK_LEN = 1024;       // length of factorial table
static unsigned long mt[N]; /* the array for the state vector  */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */

/* initializes mt[N] with a seed */
void init_genrand(unsigned long s)
{
    mt[0]= s & 0xffffffffUL;
    for (mti=1; mti<N; mti++) {
	mt[mti] =
	    (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
	/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
	/* In the previous versions, MSBs of the seed affect   */
	/* only MSBs of the array mt[].                        */
	/* 2002/01/09 modified by Makoto Matsumoto             */
	mt[mti] &= 0xffffffffUL;
	/* for >32 bit machines */
    }
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
/* slight change for C++, 2004/2/26 */
void init_by_array(unsigned long init_key[], int key_length)
{
    int i, j, k;
    init_genrand(19650218UL);
    i=1; j=0;
    k = (N>key_length ? N : key_length);
    for (; k; k--) {
	mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
	  + init_key[j] + j; /* non linear */
	mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
	i++; j++;
	if (i>=N) { mt[0] = mt[N-1]; i=1; }
	if (j>=key_length) j=0;
    }
    for (k=N-1; k; k--) {
	mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
	  - i; /* non linear */
	mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
	i++;
	if (i>=N) { mt[0] = mt[N-1]; i=1; }
    }

    mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
}

/* generates a random number on [0,0xffffffff]-interval */
unsigned long genrand_int32(void)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */

    if (mti >= N) { /* generate N words at one time */
	int kk;
	 srand((unsigned) time(NULL));

	if (mti == N+1)   /* if init_genrand() has not been called, */
	    init_genrand(rand()); // !!!!!!/* random initial seed is used */

	for (kk=0;kk<N-M;kk++) {
	    y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
	    mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
	}
	for (;kk<N-1;kk++) {
	    y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
	    mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
	}
	y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
	mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

	mti = 0;
    }

    y = mt[mti++];

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
}

/* generates a random number on [0,0x7fffffff]-interval */
long genrand_int31(void)
{
    return (long)(genrand_int32()>>1);
}

/* generates a random number on [0,1]-real-interval */
double genrand_real1(void)
{
    return genrand_int32()*(1.0/4294967295.0);
    /* divided by 2^32-1 */
}

/* generates a random number on [0,1)-real-interval */
double genrand_real2(void)
{
    return genrand_int32()*(1.0/4294967296.0);
    /* divided by 2^32 */
}

/* generates a random number on (0,1)-real-interval */
double genrand_real3(void)
{
    return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
    /* divided by 2^32 */
}

/* generates a random number on [0,1) with 53-bit resolution*/
double genrand_res53(void)
{
    unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
    return(a*67108864.0+b)*(1.0/9007199254740992.0);
}

double LnFac(int n) {
  // log factorial function. gives natural logarithm of n!

  // define constants
  static const double        // coefficients in Stirling approximation     
    C0 =  0.918938533204672722,   // ln(sqrt(2*pi))
    C1 =  1./12., 
    C3 = -1./360.;
    // C5 =  1./1260.,  // use r^5 term if FAK_LEN < 50
    // C7 = -1./1680.;  // use r^7 term if FAK_LEN < 20
  // static variables
  static double fac_table[FAK_LEN]; // table of ln(n!):
  static int initialized = 0;   // remember if fac_table has been initialized

  if (n < FAK_LEN) {
    if (n <= 1) {
      if (n < 0) FatalError("Parameter negative in LnFac function");  
      return 0;}
    if (!initialized) { // first time. Must initialize table
      // make table of ln(n!)
      double sum = fac_table[0] = 0.;
      for (int i=1; i<FAK_LEN; i++) {
        sum += log(i);
        fac_table[i] = sum;}
      initialized = 1;}
    return fac_table[n];}
    
  // not found in table. use Stirling approximation
  double  n1, r;
  n1 = n;  r  = 1. / n1;
  return (n1 + 0.5)*log(n1) - n1 + C0 + r*(C1 + r*r*C3);}

int PoissonLow (double L) {
/*
   This subfunction generates a random variate with the poisson 
   distribution for extremely low values of L.

   The method is a simple calculation of the probabilities of x = 1
   and x = 2. Higher values are ignored.

   The reason for using this method is to avoid the numerical inaccuracies 
   in other methods.
*/   
  double d, r;
  d = sqrt(L);
  if (genrand_real1() >= d) return 0;
  r = genrand_real1() * d;
  if (r > L * (1.-L)) return 0;
  if (r > 0.5 * L*L * (1.-L)) return 1;
  return 2;}
    

int PoissonInver(double L) {
/*
   This subfunction generates a random variate with the poisson 
   distribution using inversion by the chop down method (PIN).

   Execution time grows with L. Gives overflow for L > 80.

   The value of bound must be adjusted to the maximal value of L.
*/   
  const int bound = 130;             // safety bound. Must be > L + 8*sqrt(L).
  static double p_L_last = -1.;      // previous value of L
  static double p_f0;                // value at x=0
  double r;                          // uniform random number
  double f;                          // function value
  int x;                           // return value

  if (L != p_L_last) {               // set up
    p_L_last = L;
    p_f0 = exp(-L);}                 // f(0) = probability of x=0

  while (1) {  
    r = genrand_real1();  x = 0;  f = p_f0;
    do {                        // recursive calculation: f(x) = f(x-1) * L / x
      r -= f;
      if (r <= 0) return x;
      x++;
      f *= L;
      r *= x;}                       // instead of f /= x
    while (x <= bound);}}
  
  
int PoissonRatioUniforms(double L) {
/*
   This subfunction generates a random variate with the poisson 
   distribution using the ratio-of-uniforms rejection method (PRUAt).

   Execution time does not depend on L, except that it matters whether L
   is within the range where ln(n!) is tabulated.

   Reference: E. Stadlober: "The ratio of uniforms approach for generating
   discrete random variates". Journal of Computational and Applied Mathematics,
   vol. 31, no. 1, 1990, pp. 181-189.
*/
  static double p_L_last = -1.0;            // previous L
  static double p_a;                       // hat center
  static double p_h;                       // hat width
  static double p_g;                       // ln(L)
  static double p_q;                       // value at mode
  static int p_bound;                    // upper bound
  int mode;                              // mode
  double u;                                // uniform random
  double lf;                               // ln(f(x))
  double x;                                // real sample
  int k;                                 // integer sample

  if (p_L_last != L) {
    p_L_last = L;                           // Set-up
    p_a = L + 0.5;                          // hat center
    mode = (int)L;                        // mode
    p_g  = log(L);
    p_q = mode * p_g - LnFac(mode);         // value at mode
    p_h = sqrt(SHAT1 * (L+0.5)) + SHAT2;    // hat width
    p_bound = (int)(p_a + 6.0 * p_h);}    // safety-bound

  while(1) {
    u = genrand_real1();
    if (u == 0) continue;                   // avoid division by 0
    x = p_a + p_h * (genrand_real1() - 0.5) / u;
    if (x < 0 || x >= p_bound) continue;    // reject if outside valid range
    k = (int)(x);
    lf = k * p_g - LnFac(k) - p_q;
    if (lf >= u * (4.0 - u) - 3.0) break;   // quick acceptance
    if (u * (u - lf) > 1.0) continue;       // quick rejection
    if (2.0 * log(u) <= lf) break;}         // final acceptance
  return(k);}


/***********************************************************************
		      Poisson distribution
***********************************************************************/
int Poisson (double L) {
/*
   This function generates a random variate with the poisson distribution.

   Uses inversion by chop-down method for L < 17, and ratio-of-uniforms
   method for L >= 17.

   For L < 1.E-6 numerical inaccuracy is avoided by direct calculation.
*/

  //------------------------------------------------------------------
  //                 choose method
  //------------------------------------------------------------------
  if (L < 17) {
    if (L < 1.E-6) {
      if (L == 0) return 0;
      if (L < 0) FatalError("Parameter negative in poisson function");

      //--------------------------------------------------------------
      // calculate probabilities
      //--------------------------------------------------------------
      // For extremely small L we calculate the probabilities of x = 1
      // and x = 2 (ignoring higher x). The reason for using this
      // method is to prevent numerical inaccuracies in other methods.
      //--------------------------------------------------------------
      return PoissonLow(L);}

    else {

      //--------------------------------------------------------------
      // inversion method
      //--------------------------------------------------------------
      // The computation time for this method grows with L.
      // Gives overflow for L > 80
      //--------------------------------------------------------------
      return PoissonInver(L);}}

  else {
    if (L > 2.E9) FatalError("Parameter too big in poisson function");

    //----------------------------------------------------------------
    // ratio-of-uniforms method
    //----------------------------------------------------------------
    // The computation time for this method does not depend on L.
    // Use where other methods would be slower.
    //----------------------------------------------------------------
    return PoissonRatioUniforms(L);}}


/*int main(void)
{
      clrscr();
    int i; long d,sum=0,sum1=0; double mean=0,mean2=0,disp=0;
    for(i=0;i<99;i++)
    {
    d=Poisson(6.7);
    sum+=d;
    sum1+=d*d;
      printf("%i\n",d);
     }

     mean=sum/100.0;
     mean2=sum1/100.0;
      disp=sqrt(mean2-mean);
      printf("%f\t%f\n",mean,disp);
     getch();
    return 0;
}  */

// variables used by Normal distribution
  double normal_x2;  int normal_x2_valid=0;

/***********************************************************************
		      Exponential distribution
***********************************************************************/
double ExpDist(double lambda)
{
    double ksi;
    ksi=-log(genrand_real1())/lambda;
    return ksi;
}

/***********************************************************************
		      Normal distribution
***********************************************************************/

double Normal(double m, double s) {
  // normal distribution with mean m and standard deviation s
  double normal_x1;          // first random coordinate 
  double w;                  // radius
  if (normal_x2_valid) {     // we have a valid result from last call
    normal_x2_valid = 0;
    return normal_x2 * s + m;}

  // make two normally distributed variates by Box-Muller transformation
  do {
    normal_x1 = 2. * genrand_real1() - 1.;
    normal_x2 = 2. * genrand_real1() - 1.;
    w = normal_x1*normal_x1 + normal_x2*normal_x2;}
  while (w >= 1. || w < 1E-30);
  w = sqrt(log(w)*(-2./w));
  normal_x1 *= w;  normal_x2 *= w;     // normal_x1 and normal_x2 are independent normally distributed variates
  normal_x2_valid = 1;                 // save normal_x2 for next call
  return normal_x1 * s + m;}

/***********************************************************************
                     End of program
***********************************************************************/

void EndOfProgram() {
  // This function takes care of whatever is necessary to do when the 
  // program is finished

  // It may be necessary to wait for the user to press a key
  // in order to prevent the output window from disappearing.
  // Remove the #ifdef and #endif lines to unconditionally wait for a key press;
  // Remove all three lines to not wait:
  #ifdef _GETCH_DEFINED_
  getchar();                             // wait for user to press a key
  #endif

  // It may be necessary to end the program with a linefeed:
  #if defined (__unix__) || defined (_MSC_VER)
  printf("\n");                        // end program with a linefeed
  #endif

  }


/***********************************************************************
                     Error message function
***********************************************************************/

void FatalError(char * ErrorText) {
  // This function outputs an error message and aborts the program.

  // Important: There is no universally portable way of outputting an 
  // error message. You may have to modify this function to output
  // the error message in a way that is appropriate for your system.


  // Check if FatalAppExit exists (this macro is defined in winbase.h)
  #ifdef FatalAppExit  

  // in Windows, use FatalAppExit:
  FatalAppExit(0, ErrorText);

  #else

  // in console mode, print error message  
  printf ("\n%s\n", ErrorText);
  EndOfProgram();

  #endif

  // Terminate program with error code
  exit(1);}