clear;

% --- Number of sampled points and repetitions:

M = 1000;
N = 100;

% --- Outer loop:

z = 0.0;
z2 = 0.0;

for n = 1:N;

% --- Inner loop:

m = 0;

for i = 1:M;

% --- Sample point:

x = rand;
y = rand;

% --- Test and add to count:

if (x^2 + y^2 < 1.0);
m = m + 1;
endif;

endfor;

% --- Print result:

printf("pi(%3d) = %1.5f\n", n, 4*m/M);

% --- Store result:

z = z + 4*m/M;
z2 = z2 + (4*m/M)^2;

endfor;

% --- Output:

res = z/N;
sig = sqrt((1.0/(N*(N - 1.0)))*(z2 - (1.0/N)*z*z));

printf("\nN = %d M = %d pi = %1.5f +/- %1.5f [%1.5f %1.5f]\n", N, M, res, sig, res - 1.96*sig, res + 1.96*sig);