>function map MonteCarloPI (n,plot=false) ...
$  X:=random(1,n);
$  Y:=random(1,n);
$  if plot then
$      plot2d(X,Y,>points,style=".");
$      plot2d("sqrt(1-x^2)",color=2,>add);
$  endif
$  return sum(X^2+Y^2<1)/n*4;
$endfunction
>MonteCarloPI(10^(1:7))
 [ 3.6  2.96  3.224  3.1404  3.1398  3.141548  3.1421492 ]
>pi
 3.14159265359
>MonteCarloPI(10000,true):