We definitely don’t have a shortage of methods for generating normallydistributed random numbers from a source of uniformlydistributed random numbers. One of such methods is the so called Polar Method, a variation of the BoxMuller Transform, which I already described before. You might want to take a look at it before reading this.
The Polar Method was initially proposed, it seems, by the great George Marsaglia along with T. A. Bray, in the 1964 paper A Convenient Method for Generating Normal Variables. This paper must be really good, because it is available online for the whopping price of US$35.00, or seven bucks per page! Sorry, but no thanks.
Anyway, let’s see how it works. What we’ll do here is to construct a function
that uses the Polar Method to generate normallydistributed random numbers from
a source of uniformlydistributed random numbers (namely, Math.random()
). As
in the case of the BoxMuller Transform, our RNG (Random Number Generator) will
generate two random numbers at once, but we want to return only one at a time.
Therefore, here we’ll also use two global variables to let us keep that spare
random number around until we need it:
⟨Polar Method normal RNG globals⟩ =
let polarMethodNormalRNG_haveSpare = false;
let polarMethodNormalRNG_spare = 0.0;
Our RNG has the structure shown below: it makes use of these global variables and looks like a regular RNG: takes no parameters and returns a number.
⟨Polar Method normal RNG⟩ =
⟨Polar Method normal RNG globals⟩
function polarMethodNormalRNG(): number {
⟨Polar Method normal RNG body⟩
}
On to the implementation of polarMethodNormalRNG()
: we first check if we have
that spare number available. If we do, we return it immediately.
⟨Polar Method normal RNG body⟩ =
if (polarMethodNormalRNG_haveSpare) {
polarMethodNormalRNG_haveSpare = false;
return polarMethodNormalRNG_spare;
}
Otherwise, we reached the interesting part of this explanation: the point where we’ll actually use the Polar Method!
We first create two random numbers, u1
and u2
, taken from a uniform
distribution, so that the point (u1
, u2
) lies in the unit circle. Notice
that we do so by trial and error, rejecting points that fall outside it and
trying again:
⟨Polar Method normal RNG body⟩ +=
let u0, u1, s: number;
do {
u0 = Math.random() * 2.0  1.0;
u1 = Math.random() * 2.0  1.0;
s = u0*u0 + u1*u1;
} while (s >= 1  s == 0);
Then we scale this point by a magic value m
^{1} that will leave
the resulting point (z0
, z1
) in exactly the same position as the tip of the
vector we generated in the BoxMuller Transform:
⟨Polar Method normal RNG body⟩ +=
let m = Math.sqrt(2.0 * Math.log(s) / s);
let z0 = u0 * m;
let z1 = u1 * m;
With the BoxMuller Transform, the generated vector was in polar coordinates and
therefore we had to use sine and cosine to project it into the Cartesian axes
and thus obtain the desired normallydistributed numbers. But here things are
different: the point (z0
, z1
) is already born in rectangular coordinates, so
we can just use it directly^{2}: z0
and z1
are the numbers we
want. So we return z0
while keeping z1
at hand for future invocations:
⟨Polar Method normal RNG body⟩ +=
polarMethodNormalRNG_haveSpare = true;
polarMethodNormalRNG_spare = z1;
return z0;
Again, my limited mathematical knowledge is not enough to explain all the
details about how the Polar Method works.^{3} But, really, the thing
about the Polar Method is not mathematical (the math here is pretty much the
same as in the BoxMuller transform). The important idea here is computational:
even though the Polar Method discards and regenerates more than 20% of all
generated points^{4} in that while
loop, it will still be
faster than BoxMuller because it avoids calling sin()
and cos()
. (Ahem,
usually it will be faster. Don’t miss the performance
considerations below.)
Visualization
This visualization shows the Polar Method doing its work. Depending on the button you click, it will generate one, ten, or one hundred pairs of random numbers. I suggest you to first generate one point at a time to get a feeling of what the algorithm does in a microlevel; then, generate lots of points at once to see how it behaves statistically. There are some notes after the diagram that may help you understand it.
Number of discarded points:
Number of points generated:
Ratio:
What’s happening here?
 Each click in “One” shows the steps
polarMethodNormalRNG()
goes through to generate the random numbers: First it shows the point (
u0
,u1
). If the point lies outside the unit circle, the point is drawn in red and the visualization ends. (The actual algorithm, as we’ve seen, would continue trying until getting a point that lies inside the unit circle.)  If the point is inside the unit circle, it is drawn in blue and translated
to its final position: (
z0
,z1
). This movement depicts the multiplication bym
.z0
andz1
are the two random numbers generated bypolarMethodNormalRNG()
, so the Polar Method is complete at this point.
 First it shows the point (
 For clarity sake, we show one additional animated step: the point (
z0
,z1
) being projected into the vertical and horizontal axes, and the creation of histograms showing where these projections fall (in other words, the histograms show the distribution of values generated by the Polar Method).  After you generate enough pairs of numbers (say, some hundreds of pairs), both histograms should resemble a normal distribution.
 This shows that
z0
andz1
follow a normal distribution.  And that’s exactly what the Polar Method does: it transforms random numbers in a uniform distribution into random numbers in a normal distribution.
 You may be thinking: a normal distribution is unbounded, so the Polar Method could produce numbers from ∞ to +∞. How come in this visualization all the output values fit neatly in the [3, +3] range? Simple: because my visualization discards any value that falls outside of this range. This is just for the sake of a clearer visualization. I am not really cheating.
Performance Considerations
Descriptions of the Polar Method often praise it for its improved performance over the BoxMuller Transform. After all, given that the Polar Method avoids calling those expensive trigonometric functions it must be faster, right?
Well, I guess this was close an absolute truth when the Polar Method was introduced in the 1960s, but nowadays in addition to the CPUs we all love, we also have things like GPUs and FPGAs. In a nice 2009 paper called (hold you breath) A Comparison of CPUs, GPUs, FPGAs, and Massively Parallel Processor Arrays for Random Number Generation, the authors investigate different algorithms related with random number generation and how they perform in different categories of hardware.
An interesting conclusion of the paper is that for each of the four hardware
categories analyzed, a different algorithm was deemed the best. For GPUs, in
particular, the BoxMuller Transform was the top choice, despite all that
trigonometry. This happened for two reasons. First, GPUs can compute sines and
cosines really fast. Second, to make a long story short, GPUs don’t deal well
with branching and thus that while
loop in the Polar Method was a performance
disadvantage.
So, I guess it is fair to say that the Polar Method is faster than the BoxMuller transform on most of the hardware we use most of the time^{5} – but keep your mind open.
Nutrition Facts
AKA: Marsaglia Polar Method.
See also: This method is a variation of the BoxMuller Transform.
Keywords: Normal distribution, Gaussian distribution.
Source code: polar_method.ts
(the code
generated from the snippets on this page);
dump_polar_method.ts
(the interactive
diagram above); dump.ts
(helper utilities).
References
 David B. Thomas, Lee Howes, Wayne Luk. A Comparison of CPUs, GPUs, FPGAs, and Massively Parallel Processor Arrays for Random Number Generation. FPGA ‘09: Proceedings of the ACM/SIGDA International Symposium on Field Programmable Gate Arrays, 2009.
 Wikipedia, Marsaglia polar method.

I assume that
m
does not stand for “magic”, and that people with better mathematics knowledge can even explain what this number means. To me, it is magic. ↩︎ 
Which makes me wonder why this is called the Polar Method. Wouldn’t Rectangular Method be more descriptive?! (This is an honest doubt I have, by the way.) ↩︎

Hopefully, at least the visualization will help you to get an intuitive understanding of how the Polar Method works. ↩︎

We generate points in a square of area 2²=4 and discard points that fall outside an inscribed circle of area π×1²=π. So, the proportion of discarded points will be (4π)/4 ≈ 0.2146. ↩︎

Here I am required to point out the existence of the Ziggurat Algorithm. It is even faster than the Polar Method on regular CPUs (it won in the “CPU category” on the cited paper), but also has the same branching issue when running on GPUs. ↩︎