% page2p3.m
N=10000000
S1=random('unif',0,1,1,N);
S2=random('unif',0,1,1,N);

r0=0.5
x0=0.5
y0=0.5

for i0=1:10
    Nhit(i0)=0;
    for i=1:N/10
        nu=(i0-1)*(N/10)+i;
        if (S1(nu)-x0)^2+(S2(nu)-y0)^2 <= r0^2
            Nhit(i0)=Nhit(i0)+1;
        end
    end
    Nhit(i0)=Nhit(i0)/(N/10);
end

Nhit_a=0;
for i0=1:10
    Nhit_a = Nhit_a + Nhit(i0);
end
Nhit_a=Nhit_a/10;
'average number:',Nhit_a
Nhit_d=0;
for i0=1:10
    Nhit_d = Nhit_d + (Nhit(i0)-Nhit_a)^2;
end
Nhit_d=sqrt(Nhit_d)/9;
'with estimated deviation of: ',Nhit_d
'while the correct answer is: ',pi*0.5^2