#include <stdlib.h>
#include <math.h>
/* Uses the polar form of the Box-Muller transform to generate
* a random number with the given mean & standard deviation.
*
* Note that this uses the common yet nonstandard function random().
*
* Be sure to call srandom() beforehand!
*/
double bellrand(double mean, double stddev) {
double u, v, r;
do {
u = (random() % 1001) / 1000.0 * 2.0 - 1.0;
v = (random() % 1001) / 1000.0 * 2.0 - 1.0;
r = u*u + v*v;
} while (r == 0.0 || r >= 1.0);
return mean + stddev * u * sqrt(-2 * log(r) / r);
}