x <- runif(1e4, -1, 1)
 y <- runif(1e4, -1, 1)

 s <- x^2 + y^2
 x <- x[s < 1]
 y <- y[s < 1]
 s <- x^2 + y^2

##### s is Unif(0, 1)

 ss <- sort(s)
 pp <- ppoints(length(ss))

 hist(ss)
 plot(pp, ss, xlab = "quantiles of U(0, 1)", ylab = "quantiles of data")
 abline(0, 1)

##### t = - 2 log(s) is exp(1 / 2) = chisq(2)

 tt <- rev(- 2 * log(ss))
 plot(qchisq(pp, 2), tt, xlab = "quantiles of chisq(2)",
     ylab = "quantiles of data")
 abline(0, 1)

##### t  has distribution of length of bivariate std. normal r. v.
##### multiply sqrt(t) by unit vector in the (x, y) direction,
##### which is (x, y) / sqrt(s) to get bivariate std. normal r. v.

 u <- sqrt(- 2 * log(s) / s)
 xx <- x * u
 yy <- y * u
 plot(x, y, asp = 1)
 plot(xx, yy, asp = 1)
 plot(qnorm(pp), sort(xx), xlab = "quantiles of N(0, 1)",
     ylab = "quantiles of data")
 abline(0, 1)
 plot(qnorm(pp), sort(yy), xlab = "quantiles of N(0, 1)",
     ylab = "quantiles of data")
 abline(0, 1)