From: Luben
Subject: C-FAQ, q. 13.20
Message-ID: <Pine.GSO.3.95.980510025115.27152A-100000@chimp>
Date: Sun, 10 May 1998 03:12:21 -0400

The following version uses explictly the inverse of the density function of the Normal distribution, as well as you can specify mean and standard deviation. It's no longer than the one at q. 13.20 and simpler. Check it out.

/* Generate random numbers with mean mean, standard deviation std_dev, */
/* with a Normal (Gaussian) distribution.                              */
/* mean is any real number, std_dev > 0 (or machine epsilon).          */
/* The formula used is the inverse of the density function of the      */
/* Normal distribution:                                                */
/* x = mean +/- std_dev * sqrt((-2.0) * log(y)), 0 < y <= 1            */
/* The client should call srand(int) to initialize the seed,           */
/* before any calls to gauss_rand().                                   */
/* L.T., May 3, 1998, University of Toronto                            */

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

double gauss_rand(double mean, double std_dev)
{
  double   y;
  unsigned bin;

  errno = 0;

  /* std_dev must be greater than 0.0 (or machine epsilon) */
  if ( ! (std_dev > DBL_EPSILON)) {
    errno = EDOM;                      /* domain error */
    perror("gaussrand: std_dev too small or zero");
    return mean;
  }

  y = (double) rand() / (RAND_MAX + 1.0);   /* 0.0 <= y < 1.0 */
  bin  = (y < 0.5) ? 0 : 1;
  y = fabs(y - 1.0);                        /* 0.0 < y <= 1.0 */
  y = std_dev * sqrt((-2.0) * log(y));

  return bin ? (mean + y) : (mean - y);
}