rnd.c

/* @(#)rnd.c	2.1.8.2 
 *
 * 
 * RANDOM.C -- Implements Park & Miller's "Minimum Standard" RNG
 * 
 * (Reference:  CACM, Oct 1988, pp 1192-1201)
 * 
 * NextRand:  Computes next random integer
 * UnifInt:   Yields an long uniformly distributed between given bounds 
 * UnifReal: ields a real uniformly distributed between given bounds   
 * Exponential: Yields a real exponentially distributed with given mean
 * 
 */

#include "config.h"
#include <stdio.h>
#include <math.h>
#include "dss.h"
#include "rnd.h" 

char *env_config PROTO((char *tag, char *dflt));
void NthElement(long, long *);

void
dss_random(long *tgt, long lower, long upper, long stream)
{
	*tgt = UnifInt((long)lower, (long)upper, (long)stream);
	Seed[stream].usage += 1;

	return;
}

void
row_start(int t)	\
{
	int i;
	for (i=0; i <= MAX_STREAM; i++) 
		Seed[i].usage = 0 ; 
	
	return;
}

void
row_stop(int t)	\
	{ 
	int i;
	
	/* need to allow for handling the master and detail together */
	if (t == ORDER_LINE)
		t = ORDER;
	if (t == PART_PSUPP)
		t = PART;
	
	for (i=0; i <= MAX_STREAM; i++)
		if ((Seed[i].table == t) || (Seed[i].table == tdefs[t].child))
			{ 
			if (set_seeds && (Seed[i].usage > Seed[i].boundary))
				{
				fprintf(stderr, "\nSEED CHANGE: seed[%d].usage = %d\n", 
					i, Seed[i].usage); 
				Seed[i].boundary = Seed[i].usage;
				} 
			else 
				{
				NthElement((Seed[i].boundary - Seed[i].usage), &Seed[i].value);
				}
			} 
		return;
	}

void
dump_seeds(int tbl)
{
	int i;

	for (i=0; i <= MAX_STREAM; i++)
		if (Seed[i].table == tbl)
			printf("%d:\t%ld\n", i, Seed[i].value);
	return;
}

/******************************************************************

   NextRand:  Computes next random integer

*******************************************************************/

/*
 * long NextRand( long nSeed )
 */
long
NextRand(long nSeed)

/*
 * nSeed is the previous random number; the returned value is the 
 * next random number. The routine generates all numbers in the 
 * range 1 .. nM-1.
 */

{

    /*
     * The routine returns (nSeed * nA) mod nM, where   nA (the 
     * multiplier) is 16807, and nM (the modulus) is 
     * 2147483647 = 2^31 - 1.
     * 
     * nM is prime and nA is a primitive element of the range 1..nM-1.  
     * This * means that the map nSeed = (nSeed*nA) mod nM, starting 
     * from any nSeed in 1..nM-1, runs through all elements of 1..nM-1 
     * before repeating.  It never hits 0 or nM.
     * 
     * To compute (nSeed * nA) mod nM without overflow, use the 
     * following trick.  Write nM as nQ * nA + nR, where nQ = nM / nA 
     * and nR = nM % nA.   (For nM = 2147483647 and nA = 16807, 
     * get nQ = 127773 and nR = 2836.) Write nSeed as nU * nQ + nV, 
     * where nU = nSeed / nQ and nV = nSeed % nQ.  Then we have:
     * 
     * nM  =  nA * nQ  +  nR        nQ = nM / nA        nR < nA < nQ
     * 
     * nSeed = nU * nQ  +  nV       nU = nSeed / nQ     nV < nU
     * 
     * Since nA < nQ, we have nA*nQ < nM < nA*nQ + nA < nA*nQ + nQ, 
     * i.e., nM/nQ = nA.  This gives bounds on nU and nV as well:   
     * nM > nSeed  =>  nM/nQ * >= nSeed/nQ  =>  nA >= nU ( > nV ).
     * 
     * Using ~ to mean "congruent mod nM" this gives:
     * 
     * nA * nSeed  ~  nA * (nU*nQ + nV)
     * 
     * ~  nA*nU*nQ + nA*nV
     * 
     * ~  nU * (-nR)  +  nA*nV      (as nA*nQ ~ -nR)
     * 
     * Both products in the last sum can be computed without overflow   
     * (i.e., both have absolute value < nM) since nU*nR < nA*nQ < nM, 
     * and  nA*nV < nA*nQ < nM.  Since the two products have opposite 
     * sign, their sum lies between -(nM-1) and +(nM-1).  If 
     * non-negative, it is the answer (i.e., it's congruent to 
     * nA*nSeed and lies between 0 and nM-1). Otherwise adding nM 
     * yields a number still congruent to nA*nSeed, but now between 
     * 0 and nM-1, so that's the answer.
     */

    long            nU, nV;

    nU = nSeed / nQ;
    nV = nSeed - nQ * nU;       /* i.e., nV = nSeed % nQ */
    nSeed = nA * nV - nU * nR;
    if (nSeed < 0)
        nSeed += nM;
    return (nSeed);
}

/******************************************************************

   UnifInt:  Yields an long uniformly distributed between given bounds

*******************************************************************/

/*
 * long UnifInt( long nLow, long nHigh, long nStream )
 */
long
UnifInt(long nLow, long nHigh, long nStream)

/*
 * Returns an integer uniformly distributed between nLow and nHigh, 
 * including * the endpoints.  nStream is the random number stream.   
 * Stream 0 is used if nStream is not in the range 0..MAX_STREAM.
 */

{
    double          dRange;
    long            nTemp;

    if (nStream < 0 || nStream > MAX_STREAM)
        nStream = 0;

    if (nLow > nHigh)
    {
        nTemp = nLow;
        nLow = nHigh;
        nHigh = nTemp;
    }

    dRange = DOUBLE_CAST (nHigh - nLow + 1);
    Seed[nStream].value = NextRand(Seed[nStream].value);
    nTemp = (long) (((double) Seed[nStream].value / dM) * (dRange));
    return (nLow + nTemp);
}



/******************************************************************

   UnifReal:  Yields a real uniformly distributed between given bounds

*******************************************************************/

/*
 * double UnifReal( double dLow, double dHigh, long nStream )
 */
double
UnifReal(double dLow, double dHigh, long nStream)

/*
 * Returns a double uniformly distributed between dLow and dHigh,   
 * excluding the endpoints.  nStream is the random number stream.   
 * Stream 0 is used if nStream is not in the range 0..MAX_STREAM.
 */

{
    double          dTemp;

    if (nStream < 0 || nStream > MAX_STREAM)
        nStream = 0;
    if (dLow == dHigh)
        return (dLow);
    if (dLow > dHigh)
    {
        dTemp = dLow;
        dLow = dHigh;
        dHigh = dTemp;
    }
    Seed[nStream].value = NextRand(Seed[nStream].value);
    dTemp = ((double) Seed[nStream].value / dM) * (dHigh - dLow);
    return (dLow + dTemp);
}



/******************************************************************%

   Exponential:  Yields a real exponentially distributed with given mean

*******************************************************************/

/*
 * double Exponential( double dMean, long nStream )
 */
double
Exponential(double dMean, long nStream)

/*
 * Returns a double uniformly distributed with mean dMean.  
 * 0.0 is returned iff dMean <= 0.0. nStream is the random number 
 * stream. Stream 0 is used if nStream is not in the range 
 * 0..MAX_STREAM.
 */

{
    double          dTemp;

    if (nStream < 0 || nStream > MAX_STREAM)
        nStream = 0;
    if (dMean <= 0.0)
        return (0.0);

    Seed[nStream].value = NextRand(Seed[nStream].value);
    dTemp = (double) Seed[nStream].value / dM;        /* unif between 0..1 */
    return (-dMean * log(1.0 - dTemp));
}