In order to produce normally distributed random numbers from the system function rand() that outputs uniformly distributed random numbers, we will use the "transformation method" described in W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Numerical Recipes in C: The Art of Scientific Computing, (Cambridge University Press, Cambridge, 1988). The program segments below were adapted by Daniel Kiepfer from the Numerical Recipes book. The first listing is a "helper" function to produce a uniformly distributed random number between two limits from the system supplied rand() function. //returns a uniformly distributed random number between min and max double drandUniform (double min, double max) { return rand()*(max-min)/RAND_MAX+min; } Briefly, the idea of the transformation method is the following: If f(x) is a uniform distribution, how is f(y) distributed given that y=y(x)? The key is to examine the slope of the function y(x), that is, to look at y’(x). If y’(x) = dy/dx is large (in absolute value) for some values of x and y, then |dy| > |dx|. The probability of a number being in the range dx around x is the same for equal length segments dx. That is, for a uniform distribution, a random number has an equal probability of occurring per unit x, independent of x. Upon transformation, the probability of occurring in the range dx about x will be the same as the probability of occurring in the range dy about y. Now if |dy| > |dx|, that means the probability per unit y will be smaller. Also, of course, in general y’(x) will vary with x and so f(y) will not be constant. One can show: So, to produce a non-uniform distribution f(y), one just needs to find a function y(x) such that: . In principle this is all one needs, but the algorithm used below takes advantage of a technique called the Box-Muller method, in which it turns out that generating two random numbers with each function call takes less computational effort than generating just one. This explains the structure of the function in which alternate calls to the function return a value kept in a static variable gset. Further details can be found in the Numerical Recipes book. //returns a random number distributed normally, with a given SD double drandGaussian (double stdDev) { static int noExtra = 1; static double gset; double fac, r, v1, v2; if (noExtra) { do { v1 = drandUniform(-1.0, 1.0); v2 = drandUniform(-1.0, 1.0); r = v1*v1 + v2*v2; } while (r >= 1.0); fac = sqrt(-2.0*log(r)/r); gset = v1*fac; noExtra = 0; return v2*fac*stdDev; } else { noExtra = 1; return gset*stdDev; } }