#define PROGNAME "random_nums"
#define VERSION "0.4"
#define DATE  "30.10.2006"
#define AUTHOR "Nils Bluemer"

/* yields list of random numbers. */
/* NEW (0.3): Lorentz/Cauchy random numbers */
/* NEW (0.4): drand48 48-bit linear congruential generator (LCG) */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/time.h>

#ifndef PI
  #define PI 3.1415926535
#endif

void error(char error_text[])
/* standard error handler */
{
  fprintf(stderr,"\n %s run-time error\n", PROGNAME);
    fprintf(stderr,"--%s--\n",error_text);
    fprintf(stderr,"for general help use option -h\n");
    fprintf(stderr,"...now exiting to system...\n");
    exit(1);
}

void printhelp ()
{
 printf("**********************************************************\n");
 printf("%s: random number generation utility\n",PROGNAME);
 printf("Version: %s, %s by %s\n",VERSION,DATE,AUTHOR);
 printf("options: -n# number of random numbers to generate (default: 10)\n");
 printf("         -g normal distributed (gaussian) random numbers \n"); 
 printf("         -l Lorentz/Cauchy random numbers\n");
 printf("         -h this help\n");
}

int main (int argc, char *argv[])
{
  char c;
  long int i, count;
  double x,y,r,t;
  int gaussian, lorentzian;
  struct timeval start;

  count=10;
  gaussian=0;
  lorentzian=0;
/*   srand( time(0)); */
  gettimeofday(&start, 0);
/*   printf("# %d\n",start.tv_usec); */
  srand48(start.tv_usec);
  /* Zeitmessung: siehe http://www.willemer.de/informatik/cpp/timelib.htm */

  while (--argc > 0 && (*++argv)[0] == '-')
    while (c= *++argv[0])
	   switch (c) {
	   case 'n':
 	     sscanf(++argv[0],"%d\n",&count); 
             break;
           case 'g':
             gaussian=1;
             break;
           case 'l':
             lorentzian=1;
             break;
	   case 'h':
	     printhelp();
	     exit(0);
/*  	   default:  */
/*  	     error("No valid choice");  */
	   }
 
  if (gaussian==0)
    for (i=1;i<=count;i++){
/*       x=1.0*rand()/RAND_MAX; */
      x=drand48();
      if (lorentzian==1)
	x=tan(PI*x);
      printf("%f\n", x);
    }
  else
    for (i=1;i<=count/2;i++){
      /* Box Muller method */
/*       x=1.0*rand()/RAND_MAX; */
/*       y=1.0*rand()/RAND_MAX; */
      x=drand48();
      y=drand48();
      r=sqrt(-2.0*log(1-x));
      t=2*PI*y;
      printf("%f\n", r*cos(t));
      printf("%f\n", r*sin(t));
    }
}