function [ Z ] = box_muller( n )
%box_muller Box-Muller algorithm to sample from N(0,1).
%   This function implements the Box-Muller algorithm to sample from
%   N(0,1). It returns n independent standard-nomally distributed
%   pseudo-random numbers.

% number of pairs of normals to generate
m = ceil(n/2);

% generate angles and radii
theta = rand(m,1)*2*pi; % uniform on (0, 2 pi)
R = sqrt(-2*log(rand(m,1))); % Exp(1/2)

X = R .* cos(theta);
Y = R .* sin(theta);
Z = [X; Y];
Z = Z(1:n);

end