Understanding a Gaussian Sampler I recently learned that you can generate a Gaussian sampler from a uniform sampler. One such method is the Box-Muller Transform.
I naïvely implemented this transform in the following code:
//drand() returns a float in the range (0,1) uniformly
#define PI 3.1415926535897932
float normal_box_muller()
{
    return sqrt(-2.0*log(drand()))*cos(2.0*PI*drand());
}

This worked (as expected), but compared to another implementation of a Gaussian sampler, my function was slower.
//
float normal()
{
    float f, r, d1, d2;
    do {
        d1 = 2.0 * drand() - 1.0;
        d2 = 2.0 * drand() - 1.0;
        r = d1*d1 + d2*d2;
    } while (r >= 1.0 || r == 0);
    f = sqrt(-2.0 * log(r)/r);
    return d1 * f;
} 

I am having trouble understanding

*

*How this second function works. I see some similarities with Box-Muller, both have a square root term with negative 2 multiplying a natural logarithm.


*This may not be in theme with crossvalidated, but I also can't wrap my head around why this second function is 4 times as fast as my naïve Box-Muller implementation. My only guess is that computing the cosine is quite costly.

The second sampler

*

*generate a point $(X,Y)$ uniformly within the unit sphere by accept-reject based on a Uniform proposal over $(-1,1)^2$

*returns $$Z=\underbrace{\frac{X}{\sqrt{X^2+y^2}}}_\text{on the unit circle} \times \underbrace{\epsilon}_{\sim\mathcal E(1/2)\equiv\chi^2_2}$$
which is the first coordinate of a standard bivariate Normal vector. The two elements in $Z$ (radius $\epsilon$ and angle) are indeed independent. As pointed out by Glen_b, both samplers are using the polar coordinates representation of a bivariate Normal vector.
As for the difference in execution time, I am rather surprised at the advantage given to the second solution, esp. when considering it involves a while loop.
A time comparison in R gives the opposite ordering:
benchmark("resident"={
         x=rnorm(1e4)
       },
      "box-muller"={
        x=sqrt(-2*log(runif(1e4)))*cos(2*pi*runif(1e4))
       },
      "inverse"={
        x=qnorm(runif(1e4))
      },
      "circle"={
        d1=runif(1.33e4)
        r=d1^2+runif(1.33e4)^2
        x=d1[r<1]*sqrt(-2*log(r[r<1])/r[r<1])
        },
      replications = 100,
      columns = c("test", "replications", "elapsed",
                  "relative", "user.self", "sys.self"))         

with R rnorm doing slightly worse than the inverse cdf (due to a more careful inversion of the cdf to keep the relative precision high enough):
     test     replications elapsed relative user.self sys.self
    3 box-muller       100   0.089    1.534     0.089        0
    4     circle       100   0.118    2.034     0.118        0
    2    inverse       100   0.058    1.000     0.058        0
    1   resident       100   0.066    1.138     0.065        0

