/* Illustrate how to gather execution times for various   v. 990211 */
/* program components.                                              */

/* The CPU time collection code fragments should be rewritten as    */
/*   functions.                                                     */

#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>

int main() {

  long int seed = 56329;
  double rseed = 56329;

  int i;

  double   Randu(),                          /* Several different RNGs */
	   PMrandomInt(),
	   PMrandomReal(),
           urand(),
           x;
  void     seed_urand();

  long prog_begin_time, begin_time, end_time; /* to measure executiom times */
  int rv;                                     /* return value for OS call   */
  struct rusage *usage_data;    /* to measure execution time (in microsecs) */

#define NumTests 1000000

#define UserTime 1000000*(*usage_data).ru_utime.tv_sec+(*usage_data).ru_utime.tv_usec
#define SystemTime 1000000*(*usage_data).ru_stime.tv_sec+(*usage_data).ru_stime.tv_usec

/* Allocate memory for resource useage structure */
  usage_data = (struct rusage *)malloc( sizeof( struct rusage ) );

/* Record runtime to this point */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  prog_begin_time = UserTime + SystemTime;
  begin_time = prog_begin_time;

/* Generate RNs */
  for ( i = 0; i <= NumTests; i++ )
    x = Randu( &seed );

/* Report runtime for Randu */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  end_time = UserTime + SystemTime;
  printf( "  RNG    execution time (milliseconds)\n" );
  printf( "Randu:        %8d\n", (int)((end_time - begin_time)/1000.0) );

/* Record runtime to this point */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  begin_time = UserTime + SystemTime;

/* Genrate RNs */
  for ( i = 0; i <= NumTests; i++ )
    x = PMrandomInt( &seed );

/* Report runtime for PMrandomInt */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  end_time = UserTime + SystemTime;
  printf( "PMrandomInt:  %8d\n", (int)((end_time - begin_time)/1000.0) );

/* Record runtime to this point */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  begin_time = UserTime + SystemTime;

/* Genrate RNs */
  for ( i = 0; i <= NumTests; i++ )
    x = PMrandomReal( &rseed );

/* Report runtime for PMrandomReal */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  end_time = UserTime + SystemTime;
  printf( "PMrandomReal: %8d\n", (int)((end_time - begin_time)/1000.0) );

seed_urand( 56329 );

/* Record runtime to this point */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  begin_time = UserTime + SystemTime;

/* Genrate and tabulate RNs */
  for ( i = 0; i <= NumTests; i++ )
    x = urand();

/* Report runtime for urand */
  rv = getrusage( RUSAGE_SELF, usage_data );
  if ( rv == -1 )
    printf( "getrusage call failed.\n" );
  end_time = UserTime + SystemTime;
  printf( "urand:        %8d\n", (int)((end_time - begin_time)/1000.0) );
  printf( "----------------------\n" );
  printf( "Program:      %8d\n", (int)((end_time - prog_begin_time)/1000.0) );
  }

/* Bad Example!! Do not use!!                                */
/* Returns values between 0 and 1.                           */
/* Depends on 32 bit representation for integers.            */
double Randu( seed )
  long int *seed;
{
  long int a = 16807,
           mod = 2147483647;     /* 2^31 - 1 */
  double   dmod = 2147483647;    /* 2^31 - 1 */

  *seed = *seed * a;
  if ( *seed < 0 )
    *seed += mod;
  return( ( double ) *seed / dmod );
}

/* Park & Miller function (see CACM Oct. 88) integer version */
/* of their generator.                                       */
/* This should work correctly if largest integer is 2^31 - 1 */
/* or larger.                                                */
double PMrandomInt( seed )
  long int   *seed;
{
	long int  a = 16807,            /* 7^5 */
	          m = 2147483647,       /* 2^31 - 1 */
		  q = 127773,           /* m / a (int divide) */
		  r = 2836,		/* m % a */
		  lo, hi, test;
	double dm = 2147483647;

	hi = *seed / q;
	lo = *seed % q;
	test = a * lo - r * hi;
	*seed = test > 0 ? test : test + m;
	return( (double ) *seed / dm );
}

/* Park & Miller function (see CACM Oct. 88) integer version */
/* of their generator.                                       */
/* This should work correctly if largest integer is 2^31 - 1 */
/* or larger.                                                */
double PMrandomReal( seed )
  double   *seed;
{
	double  a = 16807.0,            /* 7^5 */
	        m = 2147483647.0,       /* 2^31 - 1 */
	        q = 127773.0,           /* m / a (int divide) */
	      	r = 2836.0,		/* m % a */
		lo, hi, test;
        int     inthi;

	inthi = *seed / q;
	hi = inthi;
	lo = *seed - q * hi;
	test = a * lo - r * hi;
	*seed = test > 0.0 ? test : test + m;
	return( (double ) *seed / m );
}

/* Posting from comp.sim Jan. 16, 1996          */
 
/* In 1994 George Marsaglia posted to some news groups the idea for a new
   type of random number generator he called "multiply-with-carry".  I
   have written some C code to do this, and would like to use it, but I 
   wonder if anyone out there has any experience with MWC RNGs.  Marsaglia
   has not published anything on them yet, so far as I know.
 
   They pass the usual tests, according to Marsaglia.

   Thanks very much for any info.

   Bill Simpson
 
   Code:
*/

static unsigned long int next1, next2;

static double max = 0.0;

void seed_urand( unsigned int seed1 )
 {
 next1 = seed1;
 next2 = 69069 * seed1 + 1;          /*old Marsaglia LCG*/
 return;
 }

double urand( void )
 {
 double mod =  4294967295.0; /* 2**32 - 1 */

 next1 = 30903 * ( next1 & 65535L ) + ( next1 >> 16 );
 next2 = 18000 * ( next2 & 65535L ) + ( next2 >> 16 );
 return(  ( ( next1<<16 ) + next2 ) / mod );
 }
 
/*This has period 2^59 (fast but somewhat short period, fine for me)*/
/*all bits random (unlike some LCGs, say)*/