Random Numbers: Standard Normal Distribution

     So you want to generate numbers with a standard normal distribution? Here's the code to do it:

#include <math.h>

double stdNormRand()
{
     static bool cachedAnswer = false;     // We create two at a time - we'll just remember the second for the next call
     static double y1 = -1;     // Where we store the memorized answer for every other call

     double x1, x2, w, y2;
     if (cachedAnswer)     // Check if we already memorized an answer
     {
          cachedAnswer = false;     // If we did, then next time around we need to compute 2 more
          return y1;     // Return the memorized answer!
     }
     do {
          x1 = 2.0 * (double)rand() / (double)RAND_MAX - 1.0;
          x2 = 2.0 * (double)rand() / (double)RAND_MAX - 1.0;
          w = x1 * x1 + x2 * x2;
     } while (w >= 1.0);
     w = (double) sqrt((-2.0 * log(w))/w);
     y1 = x1 * w;
     y2 = x2 * w;

     cachedAnswer = true;     // y1 now holds the next number to be returned
     return y2;
}

     So how does it work? I'm not really sure. Is it the most efficient method? I don't know. As with most of the other projects I've put on my site, it's probably not the best way to solve the problem: but it's simple to use and it works. How accurate of an approximation is it? I'm not sure - but the CDF (computed using 10,000 points) visually matched the one on wikipedia. One thing to note: it computes two at a time and remembers the second to return later. If you're doing tests where you try to exactly reproduce results by seeding the random number generator you may need to remove the memory of this function (or it could return the memorized answer from before you reset the random number generator. To do that, just remove the second to last line (chachedAnswer = true;).