Why does this algorithm generate a standard normal distribution? I have this algorithm which I encountered:
(1) Generate $U_1$, $U_2$ independently from Uniform(0,1)
(2) Set $Y_1 = -\log{U_1}, Y_2 = -\log{U_2}$. If $Y_2 > \frac{(1-Y_1)^2}{2}$, accept $(Y_1, Y_2)$. Else reject and return to step 1.
(3) Generate $U_3$ from Uniform(0,1). If $U_3 < 0.5$, accept $X=Y_1$. Else, $X=-Y_1$

Why does X follow a standard normal distribution in the end? I know $Y_1$ and $Y_2$ here are exponential R.Vs. I probably need to understand how comparing $Y_2$ and $\frac{(1-Y_1)^2}{2}$ here does the trick and the third step is probably the result of a distribution formed after step 2 which is a folded standard normal distribution.
Edit: Source: Class Notes on Monte Carlo Methods

To draw a comparison between this Normal generator (that I will consider as von Neumann's) and the Box-Müller polar generator,
#Box-Müller
bm=function(N){
  a=sqrt(-2*log(runif(N/2)))
  b=2*pi*runif(N/2)
  return(c(a*sin(b),a*cos(b)))
}

#vonNeumann
vn=function(N){
  u=-log(runif(2.64*N))
  v=-2*log(runif(2.64*N))>(u-1)^2
  w=(runif(2.64*N)<.5)-2
  return((w*u)[v])
}

here are the relative computing times
> system.time(bm(1e8))
utilisateur     système      écoulé 
     7.015       0.649       7.674 
> system.time(vn(1e8))
utilisateur     système      écoulé 
     42.483       5.713      48.222 

