normalvariate.c

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

#define NV_MAGICCONST 1.7155277699214135 /* (4 * exp(-0.5) / sqrt(2.0)); */
#define MAX_RANDOM 2147483647 /* (2 ** 31) - 1; */

/*
 * normalvariate ported from python's stdlib (random.normalvariate).
 * released under the same licence as that (python license).
 * 
 * remember to initialize the random engine before calling this function.
 */
double normalvariate(double mu, double sigma) {
    double u1, u2, z, zz;
	for (;;) {
		u1 = ((float)random()) / MAX_RANDOM;
		u2 = 1.0 - (((float)random()) / MAX_RANDOM);
		
		z = NV_MAGICCONST * (u1 - 0.5) / u2;
		zz = z * z / 4.0;
		
		if (zz <= -(log(u2))) {
			break;
		}
	}
	
	return mu + z * sigma;
}