/*
 * CS 70, Spring 2001, Assignment 8
 *
 * Implementation of random-number-generation support functions.
 *
 * This file contains two functions: setRandomSeed, which initializes
 * the drand48 seed from the system time, and expRand, which returns
 * an exponentially-distributed random number with a specified mean.
 */

#include "randomgeneration.hh"
#include <math.h>
#include <stdlib.h>
#include <sys/time.h>

void			setRandomSeed();
					// Set a "random" random seed
double			expRand(double mean);
					// Generate exponentially-distributed
					// random number

/*
 * The random-seed function needs to know the size of the machine
 * word.  It also needs to run the pseudorandom-number generator
 * (PRNG) for a few cycles to make sure it has begun a unique
 * sequence.  These two constants define those values.
 */
static const unsigned int
			LONG_BITS = sizeof(long) * 8;
					// Number of bits in a long
static const unsigned int
			RAN_ITERS = 10;	// No. of iterations for random init

/*
 * Set up a "random" random seed.  This routine works by using the
 * time of day to initialize the pseudorandom number generator.  Then,
 * since the time of day tends to be a fairly predictable number, it
 * cycles the PRNG a few times to make sure that it has entered a
 * unique part of its sequence.
 *
 * On a machine that supports /dev/random, a better way to initialize
 * the PRNG would be to read from that device.  However, doing so would
 * complicate the code without adding instructional value.
 */
void setRandomSeed()
    {
    /*
     * Get the time of day in seconds and microseconds.  Note that the
     * system clock may not have microsecond resolution, so there is
     * no guarantee that all of the microsecond bits are significant
     */
    struct timeval tv;			// Time of day, for getting usecs
    struct timezone tz;			// Dummy for gettimeofday
    (void) gettimeofday (&tv, &tz);

    /*
     * Calculate a 48-bit random seed, and use it (seed48) to
     * initialize the random-number generator.  The seed is composed
     * of the Unix time, in microsecond resolution, converted to a
     * 48-bit number.  We would have to do some pretty hairy integer
     * arithmetic with 32-bit numbers to get a correct seed value, and
     * frankly it's just not worth it.  (This is a case where access
     * to assembly language would greatly simplify the task, but that
     * would be very non-portable.)  That leaves us with two choices:
     * (1) use the "long long" type that g++ provides, but which is
     * not part of the C++ standard (yet), or (2) use double-precision
     * arithmetic to calculate intermediate results.  The latter
     * option has the advantage that a double supports 56 bits in its
     * mantissa, so it can handle the problem, but it is also
     * portable.  We'll do it that way.
     *
     * However, in using double-precision arithmetic we must be
     * careful to avoid variations in the way different machines
     * convert numbers that are larger than the maximum integer.  Some
     * truncate; others set the result to the largest integer.  So
     * we're careful to be sure that our results are never larger than
     * the largest representable integer.
     *
     * It's worth noting that the microsecond-adjusted time is
     * actually about 52 bits wide (the Unix time is up to 32 bits,
     * and the microseconds take just a hair under 20 bits).  Thus,
     * the first truncation must discard the top 4 bits.  That's not a
     * problem, since those bits are the most slowly varying of all.
     *
     * Also note that in setting seed16[2], we couldn't use a single
     * integral divisor because its value would overflow the size of a
     * long.  That's why we use a const double and divide twice.  An
     * alternative would be to type in the value of 2^32 by hand as a
     * double, which would be easy but slightly less readable.  Since
     * this routine only executes once, the time lost in doing the
     * calculation dynamically is irrelevant.
     */
    double wideSeed = (unsigned long)tv.tv_sec * 1000000.0 + tv.tv_usec;
    unsigned short seed16[3];
    const double divisor = 1ul << 16;
    unsigned long topBits = (unsigned long)(wideSeed / divisor / divisor);
    seed16[2] = (unsigned short)topBits;
    wideSeed -= topBits * divisor * divisor;
    seed16[1] = (unsigned short)(wideSeed / divisor);
    wideSeed -= seed16[1] * divisor;
    seed16[0] = (unsigned short)wideSeed;

    seed48(seed16);

    /*
     * The above calculation will tend to generate a random seed that
     * has predictable values in some bits and unpredictable values in
     * others.  Cycle the PRNG a few times so that we will have
     * departed from the predictable part of the sequence.
     */
    for (unsigned int i = 0;  i < RAN_ITERS;  i++)
	drand48();
    }

/*
 * Generate an exponentially distributed random number with a given
 * mean.  The formula used here is a well-known way of generating an
 * exponential distribution.  For more information, see Averill M. Law
 * and W. David Kelton, Simulation Modeling and Analysis, McGraw Hill,
 * 1991.  (There is a more recent 1999 edition, but I haven't yet
 * checked to be sure it contains the same material.)
 */
double expRand(
    double		mean)		// Mean of the exponential distribution
    {
    return -mean * log(drand48());
    }