%% versao 1
clear

N = 10000;
n = 0;

for i = 1 : N
    x = rand;
    y = rand;
    d = x^2 + y^2;
    if d<1
        n = n + 1;
    end
end
PIexp = 4 * n / N

%% versao 2
clear

N = 10000;
x = rand(N,1);
y = rand(N,1);
d = x .^ 2 + y .^ 2;
n = sum( d<1 );
PIexp = 4 * n / N

%% versao 3
clear
N = 1e4;
M = 100;
for i = 1 : M
    PIexp(i) = 4 * mean( sum( rand(2,N).^2 )<1 );
end
std( PIexp )

%% versao 4
clear
N = 1e4;
M = 100;
PIexp = 4 * mean( sum( rand(2,N,M).^2 )<1 );
