N = 1000

%% a typical C-like code

n_circ = 0;
for i = 1 : N
    x = rand(1);
    y = rand(1);
    D = (x.^2 + y.^2);
    if D<=1
        n_circ = n_circ+1;
    end
end
pi_mc = 4 * n_circ / N


%% matlab code - version 1

x = rand(N,1);
y = rand(N,1);
D = (x.^2 + y.^2);
pi_mc = 4 * sum( D<=1 ) / N

%% matlab code - version 2 - avoiding compute D separately

x = rand(N,1);
y = rand(N,1);
pi_mc = 4 * sum( ( x.^2 + y.^2 )<=1 ) / N

%% matlab code - version 3 - combining x and y in a two-column matrix

xy = rand(N,2);
pi_mc = 4 * sum( ( xy(:,1).^2 + xy(:,2).^2 )<=1 ) / N

%% matlab code - version 4 - combining x and y in a two-column matrix and 
%      replacing the "+" by a matricial product

xy = rand(N,2);
pi_mc = 4 * sum( ( (xy.^2)*[1; 1] )<=1 ) / N

%% matlab code - version 5 - avoiding compute x and y separately 

pi_mc = 4 * sum( ( rand(N,1).^2 + rand(N,1).^2 )<=1 ) / N


%% matlab code - version 6 - avoiding compute x and y separately and 
%     replacing the "+" by a matricial product

pi_mc = 4 * sum( ( (rand(N,2).^2)*[1; 1] )<=1 ) / N

%% matlab code - version 7 - avoiding compute x and y separately and 
%    replacing the "+" by the function "sum" over the columns (2nd dimention)

pi_mc = 4 * sum( sum( rand(N,2).^2, 2 )<=1 ) / N

%% matlab code - version 8 - avoiding compute x and y separately, replacing
%     the "+" by the function "sum" and using the function "mean" to compute
%     the fraction of the inside points

pi_mc = 4 * mean( sum( rand(N,2).^2, 2 )<=1 )

