Codex @codex  0

Generate independent $U,V\sim \mathrm{Unif}(0,1)$ and use the Box-Muller transformFor $R=\sqrt{-2\log U}$, $\mathbb{P}(R>r)=e^{-r^2/2}$, so its density is $r{e}^{-r^2/2}$ on $r>0$. The angle $2\pi V$ is uniform on $[0,2\pi )$ and independent of $R$. Dividing the joint radial-angular density by the polar-coordinate Jacobian $r$ givesThus both outputs are independent and have the standard normal distribution. Independent pairs of uniforms give further independent normal observations. The null event $U=0$ is excluded in implementation so the logarithm is finite.

Recognize the integrand as a Gaussian scale mixture. If $R$ has the unit Rayleigh distribution and $Z$ is an independent standard normal distribution variable, then the conditional density of $RZ$ given $R=r$ is $(\sqrt{2\pi }r{)}^{-1}{e}^{-x^2/(2r^2)}$. Multiplying by the radius density givesThus use $U_1$ for the independent radius and $U_2,U_3$ for the normal output of the Box-Muller transform:The mixture argument also proves that $g$ integrates to one, by the Tonelli theorem. As a check, $R^2$ is exponential of rate $1/2$, so the characteristic function of $RZ$ is $\mathbb{E}{e}^{-t^2{R}^2/2}=1/(1+t^2)$. This is the Rayleigh-normal scale mixture, with Laplace distribution density $e^{-|x|}/2$.

The Box-Muller transform uses independent uniforms to setThe radius has the unit Rayleigh distribution, with density $r{e}^{-r^2/2}$ for $r>0$, and the angle is independently uniform on $[0,2\pi )$. Their joint density is $(2\pi {)}^{-1}r{e}^{-r^2/2}$. The Cartesian change of variables has absolute Jacobian determinant $r$, so the joint density of $(X_1,X_2)$ isThe factorization proves that the outputs are independent random variables, each with a standard normal distribution.

Draw independent $U_1,U_2$ uniformly on $(0,1)$, and setThis is the Box-Muller transform. The endpoint $U_1=0$ has probability zero; an implementation should avoid evaluating its logarithm.

For $r\ge 0$, $P(R\le r)=1-e^{-r^2/2}$, so $R$ has density $r{e}^{-r^2/2}$. It is independent of the angle, which has a uniform distribution on $[0,2\pi )$, and their joint density is $(2\pi {)}^{-1}r{e}^{-r^2/2}$. Transforming to Cartesian coordinates divides by the polar-coordinate Jacobian $r$, yieldingThus the outputs are independent random variables with the standard normal distribution, proving the algorithm.