% --------------------------------------------------- 
% pi_mc.m
% Monte Carlo estimate for pi
% comment: only 2 random numbers generated at a time 
%          => slow Matlab performance
% ---------------------------------------------------
clear all; close all;
N=10;
fprintf('%20s %20s %20s\n','N','approx pi','deviation');
tab = [];
tic
while N<10^6
    Nin = 0;  
    for i=1:N
     xy = rand(1,2);     % random point 
     d  = sum(xy.^2); % distance from origin
     if (d <= 1)
       Nin = Nin + 1;
     end
    end
    piappr = 4*Nin/N;
    fprintf('%20g %20.15f %20f15\n',N,piappr,pi-piappr);
    tab = [tab;[N piappr]];
    N = 2*N;
end
semilogx(tab(:,1),tab(:,2),'bo-',tab(:,1),pi*ones(size(tab(:,1),1)),'r-');
toc
   