/*
 * randsub.c  --  Interface to random()
 * Copyright (c) 2006 by James Dow Allen
 *
 *  This program can be distributed or modified subject to
 *  the terms of the GNU General Public License (version 2
 *  or, at your option, any later version) as published by
 *  the Free Software Foundation.
 */

#include	<stdlib.h>
#include	<stdio.h>
#include	<math.h>

/* Random number seed intialization */
void seed_rand(unsigned int x)
{
	if (x == 0) {
		x = getpid() + (int)time((long *)0);
		printf("Seed for randoms = %u\n", x);
	}
	srandom(x);
}

/* Return a uniform random variable from `minv' to `maxv' */
double	ran_unif(double minv, double maxv)
{
	/* Change 30 to 14 if you don't have a 32-bit random generator */
	return minv + (maxv - minv) *
		(random() & ((1L << 30) - 1) | 1) / (float)(1L << 30);
}

/* Return a Poisson random variable with mean and variance = log(m) */
int ran_poisson(double m)
{
	static double unif_save = -1;
	double	z;
	int	n;

	z = m;
	z *= unif_save < 0 ? ran_unif(0.0, 1.0) : unif_save;
	for (n = 0; z >= 1.0; n += 1)
		z *= ran_unif(0.0, 1.0);
	unif_save = z;
	return n;
}

/* Return a gaussian random variable with mean `meanv' and std-dev `sdevv' */
double ran_gaussian(double meanv, double sdevv)
{
	static iset = 1;
	static double	v1, v2, fac, r;

	if (iset ^= 1)
		return meanv + sdevv * v1 * fac;
	do {
		v1 = ran_unif(-1.0, 1.0);
		v2 = ran_unif(-1.0, 1.0);
		r = v1 * v1 + v2 * v2;
	} while (r >= 1);
	fac = sqrt(-2 * log(r) / r);
	return meanv + sdevv * v2 * fac;
}