/*
 *
 *    Copyright 2001, Cadence Design Systems, Inc. 
 *    All rights reserved. Cadence Design Systems, 
 *    Inc. products are covered by U.S. and foreign
 *    patents, issued and pending.                 
 *                                                 
 */
/*
 * FILE: noise.c (used to be rand.c)
 * DATE: January 21, 1987 Paul "hot rod" Titchener
 * OTHER RELATED FILES: 
 * NAME:
 * PURPOSE: Set of random noise generating functions
 * DESCRIPTION:
 *  The older routines (noise*) are based on a multiplicitive
 *  congruential generator.  The newer routines (gen*) are
 *  base on the Mersenne Twister which has a MUCH longer
 *  period.  Please see mt-19937-2.c for more information.
 *
 * PROBLEMS/SUGGESTIONS:
 *
 * MODIFICATION HISTORY:
 */
#include "cgs.h"
#include "cgsnoise.h"
/*
 *  Globals and defines
 */
/*
 *  Constants for the random number generator.
 *  See D. Knuth, The Art of Computer Programming, vol. 2, p 102
 */
#define RAND_MODULUS 		4294967296L			/* The modulus for the random number generator */
#define RAND_SCALE	 	2.32830643708079737543147e-10	/* One over the modulus minus one, in float format */
#define RAND_FACTOR 		1664525L			/* The multiplicative factor for the generator */
#define RAND_OFFSET 		3L				/* The multiplicative factor for the generator */

/*
 * Functions 
 */
static void sgenrand();
static double genrand();
static double gammln();

/*
 * FUNCTION: double noiseUDN(lp_seed)
 * PURPOSE:
 *
 *  Generates uniformly distributed noise
 *  ranging from 0 -> 1. seed should be
 *  set to a non-neg number prior to the first
 *  call and then not be changed.
 *
 * DESCRIPTION:
 * RETURN VALUE:
 * INCLUDE FILES NECESSARY:
 * PROBLEMS/SUGGESTIONS:
 * GLOBALS/SIDE EFFECTS:
 */
double noiseUDN(lp_seed)
unsigned long	*lp_seed;
{

    *lp_seed = *lp_seed * RAND_FACTOR + RAND_OFFSET;

    return *lp_seed * RAND_SCALE;
}



/*
 * FUNCTION: double noiseWGN(lp_seed)
 * PURPOSE:
 *
 *   Generates white Gaussian noise 	
 *   Seed should be set and not changed.
 *
 * DESCRIPTION:
 * RETURN VALUE:
 * INCLUDE FILES NECESSARY:
 * PROBLEMS/SUGGESTIONS:
 * GLOBALS/SIDE EFFECTS:
 */
double noiseWGN(lp_seed)		
long *lp_seed;
{
	double 	x, y; /* Used to map udn to wgn	*/

	x = noiseUDN(lp_seed);
	y = noiseUDN(lp_seed);
	if (y == 0.0) y = noiseUDN(lp_seed);
	
	return cos(PI2 * x) * sqrt(-log(y*y));
}



/*
 * FUNCTION: double noiseRandBit(lp_seed, d_prob_zero)
 * PURPOSE:
 *
 *   Returns a value equal to 0 or 1.
 *   seed should be set and then unchanged.
 *
 * DESCRIPTION:
 * RETURN VALUE:
 * INCLUDE FILES NECESSARY:
 * PROBLEMS/SUGGESTIONS:
 * GLOBALS/SIDE EFFECTS:
 */
double noiseRandBit(lp_seed, d_prob_zero)
long *lp_seed;
double d_prob_zero;
{
	return (noiseUDN(lp_seed) > d_prob_zero) ? 1.0 : 0.0;
}




/*
 * FUNCTION: double noiseEXP(lp_seed)
 * PURPOSE:
 *
 *	Generates an exponentially distributed Random variable.
 *
 * DESCRIPTION:
 *
 *	Returns an exponentially distributed, positive, random deviate
 *	of unit mean, using noiseUDN(lp_seed) as the source of
 *	uniform deviates.
 *	This function was retrieved from "Numerical Recipies in C"
 *	by William Press et. al.
 *
 * RETURN VALUE:
 * INCLUDE FILES NECESSARY:
 * PROBLEMS/SUGGESTIONS:
 * GLOBALS/SIDE EFFECTS:
 */
double noiseEXP(lp_seed)
long *lp_seed;
{
	double	x;
/*
 *	Returns an exponentially distributed, positive, random 
 *	deviate of unit mean, using noiseUDN(lp_seed) as the 
 *	source of uniform deviates.
 */
	x = noiseUDN(lp_seed);
	if (x == 0.0) x = noiseUDN(lp_seed);
	
	return -log(x);
}


/*
 * FUNCTION: double noisePOISS(lp_seed)
 * PURPOSE:
 *
 *	Function for Poisson distribution
 *
 * DESCRIPTION:
 *
 *	Returns as a floating point number an integer
 *	value that is a random deviate drawn from a
 *	poisson distribution of mean xm, using a source
 *	of uniform random deviates.
 *	This function was retrieved from Chapter 7 of
 *	the "Numerical Recipies in C" by Press et. al.
 *
 * RETURN VALUE:
 * INCLUDE FILES NECESSARY:
 * PROBLEMS/SUGGESTIONS:
 * GLOBALS/SIDE EFFECTS:
 */
double noisePOISS(xm,lp_seed)
double	xm;
long *lp_seed;
{
	static double	sq, alxm, g;
	static double	oldm = -1.0;
	double		em, t, y;

	if (xm < 12.0) {
		if (xm != oldm) {
			oldm = xm;
			g = exp(-xm);
		}
		em = -1;
		t = 1.0;
		do {
			em += 1.0;
			t *= noiseUDN(lp_seed);
		} while (t > g);
	}
	else {
		if (xm != oldm) {
			oldm = xm;
			sq = sqrt(2.0 * xm);
			alxm = log(xm);
			g = xm * alxm - gammln(xm + 1.0);
		}
		do {
			do {
				y = tan(PI * noiseUDN(lp_seed));
				em = sq * y + xm;
			} while (em < 0.0);
			em = floor(em);
			t = 0.9 * (1.0 + y*y) * exp(em * alxm - gammln(em + 1.0) - g);
		} while (noiseUDN(lp_seed) > t);
	}

	return em;
}
                                                     

/*
 * FUNCTION: double gammln(lp_seed)
 * PURPOSE:
 *
 *	Function needed for Poisson distribution function
 *
 * DESCRIPTION:
 *
 *	Function gammln, retrieved from Chapter 6 of the Numerical 
 *	recipies in C book.
 *	Returns the value of the natural log of the gamma function 
 *	for xx > 0.
 *
 * RETURN VALUE:
 * INCLUDE FILES NECESSARY:
 * PROBLEMS/SUGGESTIONS:
 * GLOBALS/SIDE EFFECTS:
 */
double gammln(xx)
double xx;
{
	int		j;
	double		x, tmp, ser;
	static double	cof[6] = {76.18009173,-86.50532033,24.01409822,
				    -1.231739516,0.120858003e-2,-0.536382e-5};

	x = xx - 1.0;
	tmp = x + 5.5;
	tmp -= (x + 0.5) * log(tmp);
	ser = 1.0;
	for (j = 0; j <= 5; j++) {
	      x += 1.0;
	      ser += cof[j] / x;
	}
	return -tmp + log(2.50662827465 * ser);
}

static int seedInited = 0;
static ULong seedv;
static ULong defaultMT[625];

typedef struct seeds_s {
  ULong	seed;
  char *id;
  struct seeds_s *next;
} Seeds;
static Seeds *AllSeeds = NULL;

/*---------------------------------------------------------------
 * FUNCTION:  genRandDefaultInit
 * DESCRIPTION:
 *	void genRandDefaultInit(ul_seed)
 *
 *  This initializes the default random generator
 *
 * RETURN VALUE:  
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
void genRandDefaultInit(ul_seed)
ULong ul_seed;
{
    seedv = ul_seed;
    sgenrand(seedv, defaultMT);
    seedInited = 1;
}

/*---------------------------------------------------------------
 * FUNCTION:  genRandInit
 * DESCRIPTION:
 *	ULong *genRandInit(cp_id, cp_mode, ul_seed)
 *
 *  Returns a state vector for use in later random
 *  generator uses.  The cp_id should be the hierarchical
 *  ID of the block calling this function.  cp_mode
 *  should be "system" to use the default shared
 *  generator and should be "private" to use a private
 *  generator.  If using a private generator, ul_seed
 *  gives its seed value.  The actual seed value will
 *  be (system_seed + ul_seed + 0x50000000) to make
 *  sure that the system and private generators do
 *  not overlap, even in iterated simulations.
 *
 * RETURN VALUE:  
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
ULong *genRandInit(cp_id, cp_mode, ul_seed)
char *cp_id, *cp_mode;
ULong ul_seed;
{
    ULong *mtNew;
    if (*cp_mode == 's') {  /* system generator */
	mtNew = defaultMT;
    }
    else {
	Seeds *scan = AllSeeds;
	int i;

	while (scan) {
	    if (scan->seed == ul_seed) {
		printf("Blocks %s and %s have the same seed value of %d\n", cp_id, scan->id, ul_seed);
		printf("This will lead to these sources being correlated\n");
		exit(1);
	    }
	    scan = scan->next;
	}
	scan = (Seeds *) Oalloc(VEC_MEM, sizeof(Seeds));
	scan->id = cp_id;
	scan->seed = ul_seed;
	scan->next = AllSeeds;
	AllSeeds = scan;
	
	mtNew = (ULong *) Oalloc(VEC_MEM, sizeof(ULong) * 625);