Marsaglia polar method

The underlying theory may be summarized as follows:

If u is uniformly distributed in the interval 0 ≤ u < 1, then the point (cos(2πu), sin(2πu)) is uniformly distributed on the unit circumference x² + y² = 1, and multiplying that point by an independent random variable ρ whose distribution is

${\displaystyle \Pr(\rho <a)=\int _{0}^{a}re^{-r^{2}/2}\,dr}$

will produce a point

${\displaystyle \left(\rho \cos(2\pi u),\rho \sin(2\pi u)\right)}$

whose coordinates are jointly distributed as two independent standard normal random variables.

This idea dates back to Laplace, whom Gauss credits with finding the above

${\displaystyle I=\int _{-\infty }^{\infty }e^{-x^{2}/2}\,dx}$

by taking the square root of

${\displaystyle I^{2}=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }e^{-(x^{2}+y^{2})/2}\,dx\,dy=\int _{0}^{2\pi }\int _{0}^{\infty }re^{-r^{2}/2}\,dr\,d\theta .}$

The transformation to polar coordinates makes evident that θ is uniformly distributed (constant density) from 0 to 2π, and that the radial distance r has density

${\displaystyle re^{-r^{2}/2}.\,}$

(r² has the appropriate chi square distribution.)

This method of producing a pair of independent standard normal variates by radially projecting a random point on the unit circumference to a distance given by the square root of a chi-square-2 variate is called the polar method for generating a pair of normal random variables,

A direct application of this idea,

${\displaystyle x={\sqrt {-2\ln(u_{1})}}\cos(2\pi u_{2}),\quad y={\sqrt {-2\ln(u_{1})}}\sin(2\pi u_{2})}$

is called the Box–Muller transform, in which the chi variate is usually generated as

${\displaystyle {\sqrt {-2\ln(u_{1})}},}$

but that transform requires logarithm, square root, sine and cosine functions. On some processors, the cosine and sine of the same argument can be calculated in parallel using a single instruction.[3] Notably for Intel-based machines, one can use fsincos assembler instruction or the expi instruction (available e.g. in D), to calculate complex

${\displaystyle \operatorname {expi} (z)=e^{iz}=\cos(z)+i\sin(z),\,}$

and just separate the real and imaginary parts.

Note: To explicitly calculate the complex-polar form use the following substitutions in the general form,

Let ${\displaystyle r={\sqrt {-2\ln(u_{1})}}}$ and ${\displaystyle z=2\pi u_{2}.}$ Then

${\displaystyle \ re^{iz}={\sqrt {-2\ln(u_{1})}}e^{i2\pi u_{2}}={\sqrt {-2\ln(u_{1})}}\left[\cos(2\pi u_{2})+i\sin(2\pi u_{2})\right].}$

In contrast, the polar method here removes the need to calculate a cosine and sine. Instead, by solving for a point on the unit circle, these two functions can be replaced with the x and y coordinates normalized to the ${\displaystyle {\sqrt {x^{2}+y^{2}}}}$ radius. In particular, a random point (x, y) inside the unit circle is projected onto the unit circumference by setting ${\displaystyle s=x^{2}+y^{2}}$ and forming the point

${\displaystyle \left({\frac {x}{\sqrt {s}}},{\frac {y}{\sqrt {s}}}\right),\,}$

which is a faster procedure than calculating the cosine and sine. Some researchers argue that the conditional if instruction (for rejecting a point outside of the unit circle), can make programs slower on modern processors equipped with pipelining and branch prediction.[4] Also this procedure requires about 27% more evaluations of the underlying random number generator (only ${\displaystyle \pi /4\approx 79\%}$ of generated points lie inside of unit circle).

That random point on the circumference is then radially projected the required random distance by means of

${\displaystyle {\sqrt {-2\ln(s)}},\,}$

using the same s because that s is independent of the random point on the circumference and is itself uniformly distributed from 0 to 1.

Simple implementation in Java using the mean and standard deviation:

private static double spare;
private static boolean hasSpare = false;

public static synchronized double generateGaussian(double mean, double stdDev) {
    if (hasSpare) {
        hasSpare = false;
        return spare * stdDev + mean;
    } else {
        double u, v, s;
        do {
            u = Math.random() * 2 - 1;
            v = Math.random() * 2 - 1;
            s = u * u + v * v;
        } while (s >= 1 || s == 0);
        s = Math.sqrt(-2.0 * Math.log(s) / s);
        spare = v * s;
        hasSpare = true;
        return mean + stdDev * u * s;
    }
}

A non-thread safe implementation in C++ using the mean and standard deviation:

double generateGaussian(double mean, double stdDev) {
    static double spare;
    static bool hasSpare = false;

    if (hasSpare) {
        hasSpare = false;
        return spare * stdDev + mean;
    } else {
        double u, v, s;
        do {
            u = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0;
            v = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0;
            s = u * u + v * v;
        } while (s >= 1.0 || s == 0.0);
        s = sqrt(-2.0 * log(s) / s);
        spare = v * s;
        hasSpare = true;
        return mean + stdDev * u * s;
    }
}

C++11 GNU GCC libstdc++'s implementation of std::normal_distribution uses the Marsaglia polar method, as quoted from herein.