Marsaglia polar method

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.

A simple implementation in Python:

import math
import random


def marsaglia_sample():
    while True:
        u1 = random.uniform(-1, 1)
        u2 = random.uniform(-1, 1)
        if (w := u1**2 + u2**2) < 1:
            break
    z1 = u1 * math.sqrt(-2 * math.log(w) / w)
    z2 = u2 * math.sqrt(-2 * math.log(w) / w)
    return z1, z2

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.

A simple Julia implementation:

"""
    marsagliasample(N)

Generate `2N` samples from the standard normal distribution using the Marsaglia method.
"""
function marsagliasample(N)
    z = Array{Float64}(undef,N,2);
    for i in axes(z,1)
        s = Inf;
        while s > 1
            z[i,:] .= 2rand(2) .- 1;
            s = sum(abs2.(z[i,:]))
        end
        z[i,:] .*= sqrt(-2log(s)/s);
    end
    vec(z)
end

"""
    marsagliasample(n,μ,σ)

Generate `n` samples from the normal distribution with mean `μ` and standard deviation `σ` using the Marsaglia method.
"""
function marsagliasample(n,μ,σ)
    μ .+ σ*marsagliasample(cld(n,2))[1:n];
end

The for loop can be parallelized by using the Threads.@threads macro.