make_random_state(654321)$
set_random_state(%)$

n:10^2;
withplot:1$

/* initialize counter for points that land on the quarter disc */
a:0$
if (withplot=1) then xys:[]$

for i: 1 thru n do (
    /* square distance to origin */
    x:random(1.0),
    y:random(1.0),
    d:x^2+y^2,
    if(d<1) then a:a+1,

    if (withplot=1) then xys:append(xys,[[x,y]])
)$

a;

/* our estimate for pi */
4*a/n;
float(%);

if (withplot=1) then plot2d([[discrete,xys],[parametric,cos(t),sin(t)]],
    [style,points,lines],[t,0,%pi/2],same_xy)$