/* * 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 #include #include /* 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; }