/*************************************************************************/
/* Program file name: random_exp.c                                       */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 28, 2006                                          */
/*                                                                       */
/* (1) This C program is created for the book, "An Introduction to       */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

#include <stdio.h>
#include <time.h>
#include <math.h>
int seed;

main() {
// An example of generating an array of exponential
// random numbers.

  int n = 100;
  int i, ic;
  int t1, t2, t3, t4, t5, t6;
  time_t ts, tp;
  struct tm t;
  double rane();
  extern int seed;

// Initiate the seed from the current time
  ts = time(&tp);
  t  = *gmtime(&tp);
  t1 = t.tm_sec;
  t2 = t.tm_min;
  t3 = t.tm_hour;
  t4 = t.tm_mday;
  t5 = t.tm_mon;
  t6 = t.tm_year;
  seed = t6+70*(t5+12*(t4+31*(t3+23*(t2+59*t1))));
  if ((seed%2) == 0) seed -= 1;

  for (i=0; i<n; ++i) {
   printf("The \%dth exponential number is: %lf\n", i, rane());
  }
}

double rane() {
/* Exponential random number generator from a uniform
   random number generator. */

  double ranf();
  return -log(1-ranf());
}

double ranf() {
// Function to generate a uniform random number in [0,1]
// following x(i+1)=a*x(i) mod c with a=pow(7,5) and
// c=pow(2,31)-1.  Here the seed is a global variable.

  const int a = 16807,  c = 2147483647, q = 127773, r = 2836;
  int l, h, t;
  double cd = c;
  extern int seed;

  h = seed/q;
  l = seed%q;
  t = a*l - r*h;
  if (t > 0) seed = t;
  else seed = c + t;
  return seed/cd;
}