/*************************************************************************/
/* Program file name: random_gauss.c                                     */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 29, 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 Gaussian random
// numbers.

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

// 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-1; i+=2) {
   rang(&x,&y);
   printf("The \%dth Gaussian number is: %lf\n", i, x);
   printf("The \%dth Gaussian number is: %lf\n", i+1, y);
  }
}

void rang (x, y)
// Two Gaussian random numbers generated from two uniform
// random numbers.
double *x,*y;
{
  double pi,r1,r2;
  double ranf();
  pi =  4*atan(1);
  r1 = -log(1-ranf());
  r2 =  2*pi*ranf();
  r1 =  sqrt(2*r1);
  *x = r1*cos(r2);
  *y = r1*sin(r2);
}

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;
}