Monte Carlo methods

Monte Carlo methods
You are encouraged to solve this task according to the task description, using any language you may know.

A Monte Carlo Simulation is a way of approximating the value of a function where calculating the actual value is difficult or impossible.
It uses random sampling to define constraints on the value and then makes a sort of "best guess."

A simple Monte Carlo Simulation can be used to calculate the value for $\pi$.

If you had a circle and a square where the length of a side of the square was the same as the diameter of the circle, the ratio of the area of the circle to the area of the square would be $\pi /4$.

So, if you put this circle inside the square and select many random points inside the square, the number of points inside the circle divided by the number of points inside the square and the circle would be approximately $\pi /4$.

Task

Write a function to run a simulation like this, with a variable number of random points to select.

Also, show the results of a few different sample sizes.

For software where the number $\pi$ is not built-in, we give $\pi$ as a number of digits:

            3.141592653589793238462643383280

11l

F monte_carlo_pi(n)
   V inside = 0
   L 1..n
      V x = random:()
      V y = random:()
      I x * x + y * y <= 1
         inside++
   R 4.0 * inside / n

print(monte_carlo_pi(1000000))

Output:

3.13775

360 Assembly

*        Monte Carlo methods       08/03/2017
MONTECAR CSECT
         USING  MONTECAR,R13       base register
         B      72(R15)            skip savearea
         DC     17F'0'             savearea
         STM    R14,R12,12(R13)    save previous context
         ST     R13,4(R15)         link backward
         ST     R15,8(R13)         link forward
         LR     R13,R15            set addressability
         LA     R8,1000            isamples=1000
         LA     R6,4               i=4
       DO WHILE=(C,R6,LE,=F'7')    do i=4 to 7
         MH     R8,=H'10'            isamples=isamples*10
         ZAP    HITS,=P'0'           hits=0
         LA     R7,1                 j=1
       DO WHILE=(CR,R7,LE,R8)        do j=1 to isamples
         BAL    R14,RNDPK              call random
         ZAP    X,RND                  x=rnd
         BAL    R14,RNDPK              call random
         ZAP    Y,RND                  y=rnd
         ZAP    WP,X                   x
         MP     WP,X                   x**2
         DP     WP,ONE                 ~
         ZAP    XX,WP(8)               x**2   normalized
         ZAP    WP,Y                   y
         MP     WP,Y                   y**2
         DP     WP,ONE                 ~
         ZAP    YY,WP(8)               y**2   normalized
         AP     XX,YY                  xx=x**2+y**2
       IF CP,XX,LT,ONE THEN            if x**2+y**2<1 then
         AP     HITS,=P'1'               hits=hits+1
       ENDIF    ,                      endif
         LA     R7,1(R7)               j++
       ENDDO    ,                    enddo j
         CVD    R8,PSAMPLES          psamples=isamples
         ZAP    WP,=P'4'             4
         MP     WP,ONE               ~
         MP     WP,HITS              *hits
         DP     WP,PSAMPLES          /psamples
         ZAP    MCPI,WP(8)           mcpi=4*hits/psamples
         XDECO  R6,WC                edit i
         MVC    PG+4(1),WC+11        output i
         MVC    WC,MASK              load mask
         ED     WC,PSAMPLES          edit psamples
         MVC    PG+6(8),WC+8         output psamples
         UNPK   WC,MCPI              unpack mcpi
         OI     WC+15,X'F0'          zap sign
         MVC    PG+31(1),WC+6        output mcpi
         MVC    PG+33(6),WC+7        output mcpi decimals
         XPRNT  PG,L'PG              print buffer
         LA     R6,1(R6)             i++
       ENDDO    ,                  enddo i
         L      R13,4(0,R13)       restore previous savearea pointer
         LM     R14,R12,12(R13)    restore previous context
         XR     R15,R15            rc=0
         BR     R14                exit
RNDPK    EQU    *             ---- random number generator
         ZAP    WP,RNDSEED         w=seed    
         MP     WP,RNDCNSTA        w*=cnsta
         AP     WP,RNDCNSTB        w+=cnstb
         MVC    RNDSEED,WP+8       seed=w mod 10**15
         MVC    RND,=PL8'0'        0<=rnd<1
         MVC    RND+3(5),RNDSEED+3 return rnd
         BR     R14           ---- return
PSAMPLES DS     0D,PL8             F(15,0)
RNDSEED  DC     PL8'613058151221121'  linear congruential constant
RNDCNSTA DC     PL8'944021285986747'  "
RNDCNSTB DC     PL8'852529586767995'  "
RND      DS     PL8                fixed(15,9)
ONE      DC     PL8'1.000000000'   1 fixed(15,9)
HITS     DS     PL8                fixed(15,0)
X        DS     PL8                fixed(15,9)
Y        DS     PL8                fixed(15,9)
MCPI     DS     PL8                fixed(15,9)
XX       DS     PL8                fixed(15,9)
YY       DS     PL8                fixed(15,9)
PG       DC     CL80'10**x xxxxxxxx samples give Pi=x.xxxxxx'  buffer
MASK     DC     X'40202020202020202020202020202120'  mask CL16 15num
WC       DS     PL16               character 16
WP       DS     PL16               packed decimal 16
         YREGS
         END    MONTECAR

Output:

10**4    10000 samples give Pi=3.129200
10**5   100000 samples give Pi=3.145000
10**6  1000000 samples give Pi=3.141180
10**7 10000000 samples give Pi=3.141677

Action!

Library: Action! Tool Kit

Library: Action! Real Math

INCLUDE "H6:REALMATH.ACT"

DEFINE PTR="CARD"
DEFINE REAL_SIZE="6"
BYTE ARRAY realArray(1536)

PTR FUNC RealArrayPointer(BYTE i)
  PTR p

  p=realArray+i*REAL_SIZE
RETURN (p)

PROC InitRealArray()
  REAL r2,r255,ri,div
  REAL POINTER pow
  INT i

  IntToReal(2,r2)
  IntToReal(255,r255)

  FOR i=0 TO 255
  DO
    IntToReal(i,ri)
    RealDiv(ri,r255,div)
    pow=RealArrayPointer(i)
    Power(div,r2,pow)
  OD
RETURN

PROC CalcPi(INT n REAL POINTER pi)
  BYTE x,y
  INT i,counter
  REAL tmp1,tmp2,tmp3,r1,r4
  REAL POINTER pow
  
  counter=0
  IntToReal(1,r1)
  IntToReal(4,r4)
  
  FOR i=1 TO n
  DO
    x=Rand(0)
    pow=RealArrayPointer(x)
    RealAssign(pow,tmp1)

    y=Rand(0)
    pow=RealArrayPointer(y)
    RealAssign(pow,tmp2)

    RealAdd(tmp1,tmp2,tmp3)

    IF RealGreaterOrEqual(tmp3,r1)=0 THEN
      counter==+1
    FI
  OD

  IntToReal(counter,tmp1)
  RealMult(r4,tmp1,tmp2)
  IntToReal(n,tmp3)
  RealDiv(tmp2,tmp3,pi)
RETURN

PROC Test(INT n)
  REAL pi

  PrintF("%I samples -> ",n)
  CalcPi(n,pi)
  PrintRE(pi)
RETURN

PROC Main()
  Put(125) PutE() ;clear the screen

  PrintE("Initialization of data...")
  InitRealArray()

  Test(10)
  Test(100)
  Test(1000)
  Test(10000)
RETURN

Output:

Screenshot from Atari 8-bit computer

Initialization of data...
10 samples -> 3.2
100 samples -> 3.28
1000 samples -> 3.212
10000 samples -> 3.1156

Ada

with Ada.Text_IO;                use Ada.Text_IO;
with Ada.Numerics.Float_Random;  use Ada.Numerics.Float_Random;

procedure Test_Monte_Carlo is
   Dice : Generator;
   
   function Pi (Throws : Positive) return Float is
      Inside : Natural := 0;
   begin
      for Throw in 1..Throws loop
         if Random (Dice) ** 2 + Random (Dice) ** 2 <= 1.0 then
            Inside := Inside + 1;
         end if;
      end loop;
      return 4.0 * Float (Inside) / Float (Throws);
   end Pi;
begin
   Put_Line ("     10_000:" & Float'Image (Pi (     10_000)));
   Put_Line ("    100_000:" & Float'Image (Pi (    100_000)));
   Put_Line ("  1_000_000:" & Float'Image (Pi (  1_000_000)));
   Put_Line (" 10_000_000:" & Float'Image (Pi ( 10_000_000)));
   Put_Line ("100_000_000:" & Float'Image (Pi (100_000_000)));
end Test_Monte_Carlo;

The implementation uses built-in uniformly distributed on [0,1] random numbers. Note that the accuracy of the result depends on the quality of the pseudo random generator: its circle length and correlation to the function being simulated.

Output:

     10_000: 3.13920E+00
    100_000: 3.14684E+00
  1_000_000: 3.14197E+00
 10_000_000: 3.14215E+00
100_000_000: 3.14151E+00

ALGOL 68

Works with: ALGOL 68 version Standard - no extensions to language used

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

Works with: ELLA ALGOL 68 version Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386

PROC pi = (INT throws)REAL:
BEGIN
   INT inside := 0;
   TO throws DO
      IF random ** 2 + random ** 2 <= 1 THEN
         inside +:= 1
      FI
   OD;
   4 * inside / throws
END # pi #;

print (("     10 000:",pi (     10 000),new line));
print (("    100 000:",pi (    100 000),new line));
print (("  1 000 000:",pi (  1 000 000),new line));
print ((" 10 000 000:",pi ( 10 000 000),new line));
print (("100 000 000:",pi (100 000 000),new line))

Output:

     10 000:+3.15480000000000e  +0
    100 000:+3.12948000000000e  +0
  1 000 000:+3.14169200000000e  +0
 10 000 000:+3.14142040000000e  +0
100 000 000:+3.14153276000000e  +0

Arturo

Pi: function [throws][
    inside: new 0.0
    do.times: throws [
        if 1 > hypot random 0 1.0 random 0 1.0 -> inc 'inside
    ]
    return 4 * inside / throws
]
 
loop [100 1000 10000 100000 1000000] 'n ->
    print [pad to :string n 8 "=>" Pi n]

Output:

     100 => 3.4 
    1000 => 3.112 
   10000 => 3.1392 
  100000 => 3.14368 
 1000000 => 3.14106

AutoHotkey

Search autohotkey.com: Carlo methods
Source: AutoHotkey forum by Laszlo

MsgBox % MontePi(10000)   ; 3.154400 
MsgBox % MontePi(100000)  ; 3.142040 
MsgBox % MontePi(1000000) ; 3.142096 

MontePi(n) { 
   Loop %n% { 
      Random x, -1, 1.0 
      Random y, -1, 1.0 
      p += x*x+y*y < 1 
   } 
   Return 4*p/n 
}

AWK

# --- with command line argument "throws" ---

BEGIN{ th=ARGV[1];
 for(i=0; i<th; i++) cin += (rand()^2 + rand()^2) < 1 
 printf("Pi = %8.5f\n",4*cin/th)
}

usage: awk -f pi 2300

Pi =  3.14333

BASIC

Works with: QuickBasic version 4.5

Translation of: Java

DECLARE FUNCTION getPi! (throws!)
CLS
PRINT getPi(10000)
PRINT getPi(100000)
PRINT getPi(1000000)
PRINT getPi(10000000)

FUNCTION getPi (throws)
	inCircle = 0
		FOR i = 1 TO throws
			'a square with a side of length 2 centered at 0 has
			'x and y range of -1 to 1
			randX = (RND * 2) - 1'range -1 to 1
			randY = (RND * 2) - 1'range -1 to 1
			'distance from (0,0) = sqrt((x-0)^2+(y-0)^2)
			dist = SQR(randX ^ 2 + randY ^ 2)
			IF dist < 1 THEN 'circle with diameter of 2 has radius of 1
				inCircle = inCircle + 1
			END IF
		NEXT i
	getPi = 4! * inCircle / throws
END FUNCTION

Output:

3.16
3.13648
3.142828
3.141679

BASIC256

Works with: basic256 version 1.1.4.0

# Monte Carlo Simulator
# Determine value of pi
# 21010513


tosses = 1000
in_c = 0
i = 0

for i = 1 to tosses
     x = rand
     y = rand
     x2 = x * x
     y2 = y * y
     xy = x2 + y2
     d_xy = sqr(xy)
     if d_xy <= 1 then
         in_c += 1
     endif
next i

print float(4*in_c/tosses)

Output:

Throws     Result
1000       3.208
10000      3.142
20000      3.1388
40000      3.1452

Other solution:

print " Number of throws   Ratio (Pi)      Error"

for pow = 2 to 8
	n = 10 ^ pow
	pi_ = getPi(n)
	error_ = 3.141592653589793238462643383280 - pi_
	print rjust(string(int(n)), 17); "   "; ljust(string(pi_), 13); "   "; ljust(string(error_), 13)
next
end

function getPi(n)
	incircle = 0.0
	for throws = 0 to n
		incircle = incircle + (rand()^2 + rand()^2 < 1)
	next
	return 4.0 * incircle / throws
end function

Output:

 Number of throws   Ratio (Pi)      Error
              100   2.970297        0.17129562389
             1000   3.14085914086   0.00073351273
            10000   3.13208679132   0.00950586227
           100000   3.14428855711   -0.00269590352
          1000000   3.14041685958   0.00117579401
         10000000   3.14094968591   0.000643     
        100000000   3.14153         0.00006264501

BBC BASIC

      PRINT FNmontecarlo(1000)
      PRINT FNmontecarlo(10000)
      PRINT FNmontecarlo(100000)
      PRINT FNmontecarlo(1000000)
      PRINT FNmontecarlo(10000000)
      END
      
      DEF FNmontecarlo(t%)
      LOCAL i%, n%
      FOR i% = 1 TO t%
        IF RND(1)^2 + RND(1)^2 < 1 n% += 1
      NEXT
      = 4 * n% / t%

Output:

     3.136
    3.1396
   3.13756
  3.143624
 3.1412816

FreeBASIC

' version 23-10-2016
' compile with: fbc -s console

Randomize Timer  'seed the random function

Dim As Double x, y, pi, error_
Dim As UInteger m = 10, n, n_start, n_stop = m, p

Print
Print " Mumber of throws  Ratio (Pi)     Error"
Print

Do
    For n = n_start To n_stop -1
        x = Rnd
        y = Rnd
        If (x * x + y * y) <= 1 Then p = p +1
    Next
    Print Using "    ############,  "; m ;
    pi = p * 4 / m
    error_ = 3.141592653589793238462643383280 - pi
    Print RTrim(Str(pi),"0");Tab(35); Using "##.#############"; error_
    m = m * 10
    n_start = n_stop
    n_stop = m
Loop Until m > 1000000000 ' 1,000,000,000


' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End

Output:

 Mumber of throws  Ratio (Pi)     Error

               10  3.2            -0.0584073464102
              100  3.16           -0.0184073464102
            1,000  3.048           0.0935926535898
           10,000  3.1272          0.0143926535898
          100,000  3.13672         0.0048726535898
        1,000,000  3.14148         0.0001126535898
       10,000,000  3.1417668      -0.0001741464102
      100,000,000  3.14141         0.0001826535898
    1,000,000,000  3.14169192     -0.0000992664102

Liberty BASIC

 for pow = 2 to 6
    n = 10^pow
    print n, getPi(n)
next

end

function getPi(n)
    incircle = 0
    for throws=0 to n
        scan
        incircle = incircle + (rnd(1)^2+rnd(1)^2 < 1)
    next
    getPi = 4*incircle/throws
end function

Output:

100           2.89108911
1000          3.12887113
10000         3.13928607
100000        3.13864861
1000000       3.13945686

Locomotive Basic

10 mode 1:randomize time:defint a-z
20 input "How many samples";n
30 u=n/100+1
40 r=100
50 for i=1 to n
60 if i mod u=0 then locate 1,3:print using "##% done"; i/n*100
70 x=rnd*2*r-r
80 y=rnd*2*r-r
90 if sqr(x*x+y*y)<r then m=m+1
100 next
110 pi2!=4*m/n
120 locate 1,3
130 print m;"points in circle"
140 print "Computed value of pi:"pi2!
150 print "Difference to real value of pi: ";
160 print using "+#.##%"; (pi2!-pi)/pi*100

Run BASIC

Translation of: Liberty BASIC

for pow = 2 to 6
    n = 10 ^ pow
    print n; chr$(9); getPi(n)
next 
end

function getPi(n)
    incircle = 0
    for throws = 0 to n
        incircle = incircle + (rnd(1)^2 + rnd(1)^2 < 1)
    next
    getPi = 4 * incircle / throws
end function

Output:

100	3.12
1000	3.108
10000	3.1652
100000	3.14248
1000000	3.1435

True BASIC

Translation of: BASIC

FUNCTION getpi(throws)
    LET incircle = 0
    FOR i = 1 to throws
        !a square with a side of length 2 centered at 0 has
        !x and y range of -1 to 1
        LET randx = (rnd*2)-1     !range -1 to 1
        LET randy = (rnd*2)-1     !range -1 to 1
        !distance from (0,0) = sqrt((x-0)^2+(y-0)^2)
        LET dist = sqr(randx^2+randy^2)
        IF dist < 1 then          !circle with diameter of 2 has radius of 1
           LET incircle = incircle+1
        END IF
    NEXT i
    LET getpi = 4*incircle/throws
END FUNCTION

CLEAR
PRINT getpi(10000)
PRINT getpi(100000)
PRINT getpi(1000000)
PRINT getpi(10000000)
END

Output:

3.1304
3.14324
3.141716
3.1416452

PureBasic

OpenConsole()
 
Procedure.d MonteCarloPi(throws.d)
	inCircle.d = 0
		For i = 1 To throws.d
			randX.d = (Random(2147483647)/2147483647)*2-1
			randY.d = (Random(2147483647)/2147483647)*2-1 
			dist.d  = Sqr(randX.d*randX.d + randY.d*randY.d)
			If dist.d < 1 
				inCircle = inCircle + 1
			EndIf
		Next i
	pi.d = (4 * inCircle / throws.d)	
	ProcedureReturn pi.d
	
EndProcedure

PrintN ("'built-in' #Pi         = " + StrD(#PI,20))
PrintN ("MonteCarloPi(10000)    = " + StrD(MonteCarloPi(10000),20))
PrintN ("MonteCarloPi(100000)   = " + StrD(MonteCarloPi(100000),20))
PrintN ("MonteCarloPi(1000000)  = " + StrD(MonteCarloPi(1000000),20))
PrintN ("MonteCarloPi(10000000) = " + StrD(MonteCarloPi(10000000),20))

PrintN("Press any key"): Repeat: Until Inkey() <> ""

Output:

'built-in' #PI         = 3.14159265358979310000
MonteCarloPi(10000)    = 3.17119999999999980000
MonteCarloPi(100000)   = 3.14395999999999990000
MonteCarloPi(1000000)  = 3.14349599999999980000
MonteCarloPi(10000000) = 3.14127720000000020000
Press any key

C

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
double pi(double tolerance)
{
	double x, y, val, error;
	unsigned long sampled = 0, hit = 0, i;

	do {
		/* don't check error every turn, make loop tight */
		for (i = 1000000; i; i--, sampled++) {
			x = rand() / (RAND_MAX + 1.0);
			y = rand() / (RAND_MAX + 1.0);
			if (x * x + y * y < 1) hit ++;
		}

		val = (double) hit / sampled;
		error = sqrt(val * (1 - val) / sampled) * 4;
		val *= 4;

		/* some feedback, or user gets bored */
		fprintf(stderr, "Pi = %f +/- %5.3e at %ldM samples.\r",
			val, error, sampled/1000000);
	} while (!hit || error > tolerance);
              /* !hit is for completeness's sake; if no hit after 1M samples,
                 your rand() is BROKEN */

	return val;
}

int main()
{
	printf("Pi is %f\n", pi(3e-4)); /* set to 1e-4 for some fun */
	return 0;
}

C#

using System;

class Program {
    static double MonteCarloPi(int n) {
        int inside = 0;
        Random r = new Random();

        for (int i = 0; i < n; i++) {
            if (Math.Pow(r.NextDouble(), 2)+ Math.Pow(r.NextDouble(), 2) <= 1) {
                inside++;
            }
        }

        return 4.0 * inside / n;
    }

    static void Main(string[] args) {
        int value = 1000;
        for (int n = 0; n < 5; n++) {
            value *= 10;
            Console.WriteLine("{0}:{1}", value.ToString("#,###").PadLeft(11, ' '), MonteCarloPi(value));
        }
    }
}

Output:

     10,000:3.1436
    100,000:3.14632
  1,000,000:3.139476
 10,000,000:3.1424476
100,000,000:3.1413976

C++

#include<iostream>
#include<math.h>
#include<stdlib.h>
#include<time.h>
 
using namespace std;
int main(){
    int jmax=1000; // maximum value of HIT number. (Length of output file)
    int imax=1000; // maximum value of random numbers for producing HITs.
    double x,y;    // Coordinates
    int hit;       // storage variable of number of HITs
    srand(time(0));
    for (int j=0;j<jmax;j++){
        hit=0;
        x=0; y=0;
        for(int i=0;i<imax;i++){
            x=double(rand())/double(RAND_MAX);
            y=double(rand())/double(RAND_MAX);
        if(y<=sqrt(1-pow(x,2))) hit+=1; }          //Choosing HITs according to analytic formula of circle
    cout<<""<<4*double(hit)/double(imax)<<endl; }  // Print out Pi number
}

Clojure

(defn calc-pi [iterations]
  (loop [x (rand) y (rand) in 0 total 1] 
    (if (< total iterations)
      (recur (rand) (rand) (if (<= (+ (* x x) (* y y)) 1) (inc in) in) (inc total))
      (double (* (/ in total) 4)))))

(doseq [x (take 5 (iterate #(* 10 %) 10))] (println (str (format "% 8d" x) ": " (calc-pi x))))

Output:

     100: 3.2
    1000: 3.124
   10000: 3.1376
  100000: 3.14104
 1000000: 3.141064

(defn experiment 
  [] 
  (if (<= (+ (Math/pow (rand) 2) (Math/pow (rand) 2)) 1) 1 0))

(defn pi-estimate 
  [n] 
  (* 4 (float (/ (reduce + (take n (repeatedly experiment))) n))))

(pi-estimate 10000)

Output:

     3.1347999572753906

Common Lisp

(defun approximate-pi (n)
  (/ (loop repeat n count (<= (abs (complex (random 1.0) (random 1.0))) 1.0)) n 0.25))

(dolist (n (loop repeat 5 for n = 1000 then (* n 10) collect n))
  (format t "~%~8d -> ~f" n (approximate-pi n)))

Output:

    1000 -> 3.132
   10000 -> 3.1184
  100000 -> 3.1352
 1000000 -> 3.142072
10000000 -> 3.1420677

Crystal

Translation of: Ruby

def approx_pi(throws)
  times_inside = throws.times.count {Math.hypot(rand, rand) <= 1.0}
  4.0 * times_inside / throws
end
 
[1000, 10_000, 100_000, 1_000_000, 10_000_000].each do |n| 
   puts "%8d samples: PI = %s" % [n, approx_pi(n)]
end

Output:

    1000 samples: PI = 3.1
   10000 samples: PI = 3.1428
  100000 samples: PI = 3.1454
 1000000 samples: PI = 3.141012
10000000 samples: PI = 3.141148

D

import std.stdio, std.random, std.math;

double pi(in uint nthrows) /*nothrow*/ @safe /*@nogc*/ {
    uint inside;
    foreach (immutable i; 0 .. nthrows)
        if (hypot(uniform01, uniform01) <= 1)
            inside++;
    return 4.0 * inside / nthrows;
}

void main() {
    foreach (immutable p; 1 .. 8)
        writefln("%10s: %07f", 10 ^^ p, pi(10 ^^ p));
}

Output:

        10: 3.200000
       100: 3.120000
      1000: 3.076000
     10000: 3.140400
    100000: 3.146520
   1000000: 3.140192
  10000000: 3.141476

Output:

with foreach(p;1..10)

        10: 3.200000
       100: 3.240000
      1000: 3.180000
     10000: 3.150400
    100000: 3.143080
   1000000: 3.140996
  10000000: 3.141442
 100000000: 3.141439
1000000000: 3.141559

More Functional Style

void main() {
    import std.stdio, std.random, std.math, std.algorithm, std.range;

    immutable isIn = (int) => hypot(uniform01, uniform01) <= 1;
    immutable pi = (in int n)  => 4.0 * n.iota.count!isIn / n;

    foreach (immutable p; 1 .. 8)
        writefln("%10s: %07f", 10 ^^ p, pi(10 ^^ p));
}

Output:

        10: 3.200000
       100: 3.320000
      1000: 3.128000
     10000: 3.140800
    100000: 3.128400
   1000000: 3.142836
  10000000: 3.141550

Dart

From example at Dart Official Website

import 'dart:async';
import 'dart:html';
import 'dart:math' show Random;

// We changed 5 lines of code to make this sample nicer on
// the web (so that the execution waits for animation frame, 
// the number gets updated in the DOM, and the program ends 
// after 500 iterations).

main() async {
  print('Compute π using the Monte Carlo method.');
  var output = querySelector("#output");
  await for (var estimate in computePi().take(500)) {
    print('π ≅ $estimate');
    output.text = estimate.toStringAsFixed(5);
    await window.animationFrame;
  }
}

/// Generates a stream of increasingly accurate estimates of π.
Stream<double> computePi({int batch: 100000}) async* {
  var total = 0;
  var count = 0;
  while (true) {
    var points = generateRandom().take(batch);
    var inside = points.where((p) => p.isInsideUnitCircle);
    total += batch;
    count += inside.length;
    var ratio = count / total;
    // Area of a circle is A = π⋅r², therefore π = A/r².
    // So, when given random points with x ∈ <0,1>,
    // y ∈ <0,1>, the ratio of those inside a unit circle
    // should approach π / 4. Therefore, the value of π
    // should be:
    yield ratio * 4;
  }
}

Iterable<Point> generateRandom([int seed]) sync* {
  final random = new Random(seed);
  while (true) {
    yield new Point(random.nextDouble(), random.nextDouble());
  }
}

class Point {
  final double x, y;
  const Point(this.x, this.y);
  bool get isInsideUnitCircle => x * x + y * y <= 1;
}

Output:

The script give in reality an output formatted in HTML

π ≅ 3.14139

Delphi

Works with: Delphi version 6.0

Library: SysUtils,StdCtrls

function MonteCarloPi(N: cardinal): double;
{Approximate Pi by seeing if points fall inside circle}
var I,InsideCnt: integer;
var X,Y: double;
begin
InsideCnt:=0;
for I:=1 to N do
	begin
	{Random X,Y = 0..1}
	X:=Random;
	Y:=Random;
	{See if it falls in Unit Circle}
	if X*X + Y*Y <= 1 then Inc(InsideCnt);
	end;
{Because X and Y are squared, they only fall with 1/4 of the circle}
Result:=4 * InsideCnt / N;
end;


procedure ShowOneSimulation(Memo: TMemo; N: cardinal);
var MyPi: double;
begin
MyPi:=MonteCarloPi(N);
Memo.Lines.Add(Format('Samples: %15.0n   Pi= %2.15f',[N+0.0,MyPi]));
end;


procedure ShowMonteCarloPi(Memo: TMemo);
begin
ShowOneSimulation(Memo,1000);
ShowOneSimulation(Memo,10000);
ShowOneSimulation(Memo,100000);
ShowOneSimulation(Memo,1000000);
ShowOneSimulation(Memo,10000000);
ShowOneSimulation(Memo,100000000);
end;

Output:

Samples:           1,000   Pi= 3.156000000000000
Samples:          10,000   Pi= 3.152000000000000
Samples:         100,000   Pi= 3.142920000000000
Samples:       1,000,000   Pi= 3.140864000000000
Samples:      10,000,000   Pi= 3.141990800000000
Samples:     100,000,000   Pi= 3.141426720000000

E

This computes a single quadrant of the described square and circle; the effect should be the same since the other three are symmetric.

def pi(n) {
    var inside := 0
    for _ ? (entropy.nextFloat() ** 2 + entropy.nextFloat() ** 2 < 1) in 1..n {
         inside += 1
    }
    return inside * 4 / n
}

Some sample runs:

? pi(10)
# value: 2.8

? pi(10)
# value: 2.0

? pi(100) 
# value: 2.96

? pi(10000)
# value: 3.1216

? pi(100000)
# value: 3.13088 

? pi(100000)
# value: 3.13848

EasyLang

func mc n .
   for i = 1 to n
      x = randomf
      y = randomf
      if x * x + y * y < 1
         hit += 1
      .
   .
   return 4.0 * hit / n
.
numfmt 4 0
print mc 10000
print mc 100000
print mc 1000000
print mc 10000000

Output:

3.1292
3.1464
3.1407
3.1413

EDSAC order code

Because real numbers on EDSAC were restricted to the interval [-1,1), this solution estimates pi/10 instead of pi. With 100,000 trials the program would have taken about 3.5 hours on the original EDSAC.

[Monte Carlo solution for Rosetta Code.]
[EDSAC program, Initial Orders 2.]

  [Arrange the storage]
          T45K P56F     [H parameter: library s/r P1 to print real number]
          T46K P78F     [N parameter: library s/r P7 to print integer]
          T47K P210F    [M parameter: main routine]
          T48K P114F    [& (delta) parameter: library s/r C6 (division)]
          T49K P150F    [L parameter: library subroutine R4 to read data]
          T51K P172F    [G parameter: generator for pseudo-random numbers]

[Library subroutine M3, runs at load time and is then overwritten.
 Prints header; here the header sets teleprinter to figures.]
      PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
      *!!!!TRIALS!!EST!PI#X10@&..
          PK            [after header, blank tape and PK (WWG, 1951, page 91)]

[================ Main routine ====================]
          E25K TM GK
[Variables]
    [0]   PF PF         [target count: print result when count = target]
    [2]   PF PF         [count of points]
    [4]   PF PF         [count of hits (point inside circle)]
    [6]   PF PF         [x-coordinate - 1/2]
[Constants]
          T8#Z PF T10#Z PF [clear sandwich bits in 35-bit constants]
          T8Z           [resume normal loading]
    [8]   PD PF         [35-bit constant 1]
   [10]   L1229F Y819F  [35-bit constant 2/5 (near enough)]
   [12]   IF            [1/2]
   [13]   RF            [1/4]
   [14]   #F            [figures shift]
   [15]   MF            [dot (decimal point) in figures mode]
   [16]   @F            [carriage return]
   [17]   &F            [line feed]
   [18]   !F            [space]

 [Enter with acc = 0]
   [19]   A19@ GL       [read seed for LCG into 0D]
          AD T4D        [pass seed to LCG in 4D]
   [23]   A23@ GG       [initialize LCG]
          T2#@ T4#@     [zero trials and hits]
[Outer loop: round target counts]
   [27]   TF            [clear acc]
   [28]   A28@ GL       [read next target count into 0D]
          SD            [acc := -target]
          E85@          [exit if target = 0]
          T#@           [store negated target]
[Inner loop : round points in the square]
   [33]   TF T4D        [pass LCG range = 0 to return random real in [0,1)]
   [35]   A35@ G1G      [call LCG, 0D := random x]
          AD S12@ T6#@  [store x - 1/2 over next call]
          T4D
   [41]   A41@ G1G      [call LCG, 0D := random y]
          AD S12@ TD    [store y - 1/2]
          H6#@ V6#@     [acc := (x - 1/2)^2]
          HD VD         [acc := acc := (x - 1/2)^2 + (y - 1/2)^2]
          S13@          [test for point inside circle, i.e. acc < 1/4]
          E56@          [skip if not]
          TF A4#@ A8#@ T4#@ [inc number of hits]
   [56]   TF A2#@ A8#@ U2#@ [inc number of trials]
          A#@           [add negated target]
          G33@          [if not reached target, loop back]
          A2#@ TD       [pass number of trials to print s/r]
   [64]   A64@ GN       [print number of trials]
          A4#@ TD A2#@ T4D [pass hits and trials to division s/r]
   [70]   A70@ G&       [0D := hits/trials, estimated value of pi/4]
          HD V10#@ TD   [times 2/5; pass estimated pi/10 to print s/r]
          O18@ O18@ O8@ O15@ [print '  0.']
   [79]   A79@ GH P5F   [print estimated pi/10 to 5 decimals]
          O16@ O17@     [print CR, LF]
          E27@          [loop back for new target]
   [85]   O14@          [exit: print dummy character to flush printer buffer]
          ZF            [halt program]

[==================== Generator for pseudo-random numbers ===========]
[Linear congruential generator, same algorithm as Delphi 7 LCG.
 38 locations]
          E25K TG
          GK G10@ G15@ T2#Z PF T2Z I514D P257F T4#Z PF T4Z PD PF T6#Z PF T6Z PF RF A6#@ S4#@ T6#@ E25F E8Z PF T8Z PF PF A3F T14@ A4D T8#@ ZF A3F T37@ H2#@ V8#@ L512F L512F L1024F A4#@ T8#@ H6#@ C8#@ T8#@ S4D G32@ TD A8#@ E35@ H4D TD V8#@ L1F TD ZF

[==================== LIBRARY SUBROUTINES ============================]
[D6: Division, accurate, fast.
 36 locations, workspace 6D and 8D.
 0D := 0D/4D, where 4D <> 0, -1.]
          E25K T& GK
          GKA3FT34@S4DE13@T4DSDTDE2@T4DADLDTDA4DLDE8@RDU4DLDA35@
          T6DE25@U8DN8DA6DT6DH6DS6DN4DA4DYFG21@SDVDTDEFW1526D

[R4: Input of one signed integer at runtime.
 22 storage locations; working positions 4, 5, and 6.]
          E25K TL
          GKA3FT21@T4DH6@E11@P5DJFT6FVDL4FA4DTDI4FA4FS5@G7@S5@G20@SDTDT6FEF

[P1: Prints non-negative fraction in 0D, without '0.']
          E25K TH
          GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F

[P7, prints long strictly positive integer;
 10 characters, right justified, padded left with spaces.
 Even address; 35 storage locations; working position 4D.]
          E25K TN
          GKA3FT26@H28#@NDYFLDT4DS27@TFH8@S8@T1FV4DAFG31@SFLDUFOFFFSF
          L4FT4DA1FA27@G11@XFT28#ZPFT27ZP1024FP610D@524D!FO30@SFL8FE22@

[===================================================================]
[The following, without the comments and white space, might have
 been input from a separate tape.]
          E25K TM GK
          E19Z          [define entry point]
          PF            [acc = 0 on entry]
[Integers supplied by user: (1) seed for LCG; (2) list of numbers of trials
 for which to print result; increasing order, terminated by 0.
 To be read by library subroutine R4; sign comes after value.]
          987654321+100+1000+10000+100000+0+

Output:

    TRIALS  EST PI/10
       100  0.32400
      1000  0.31319
     10000  0.31371
    100000  0.31410

Elixir

defmodule MonteCarlo do
  def pi(n) do
    count = Enum.count(1..n, fn _ ->
      x = :rand.uniform
      y = :rand.uniform
      :math.sqrt(x*x + y*y) <= 1
    end)
    4 * count / n
  end
end

Enum.each([1000, 10000, 100000, 1000000, 10000000], fn n ->
  :io.format "~8w samples: PI = ~f~n", [n, MonteCarlo.pi(n)]
end)

Output:

    1000 samples: PI = 3.112000
   10000 samples: PI = 3.127200
  100000 samples: PI = 3.145440
 1000000 samples: PI = 3.142904
10000000 samples: PI = 3.141124

Erlang

With inline test

-module(monte).
-export([main/1]).

monte(N)->
    monte(N,0,0).

monte(0,InCircle,NumPoints) ->
    4 * InCircle / NumPoints;

monte(N,InCircle,NumPoints)->
    Xcoord = rand:uniform(),
    Ycoord = rand:uniform(),
    monte(N-1,
          if Xcoord*Xcoord + Ycoord*Ycoord < 1 -> InCircle + 1; true -> InCircle end, 
          NumPoints + 1).

main(N) -> io:format("PI: ~w~n", [ monte(N) ]).

Output:

8> [monte:main(X) || X <- [10000,100000,100000,10000000] ].     
PI: 3.136
PI: 3.1464
PI: 3.1412
PI: 3.1416704
[ok,ok,ok,ok]

With test in a function

-module(monte2).
-export([main/1]).

monte(N)->
    monte(N,0,0).

monte(0,InCircle,NumPoints) ->
    4 * InCircle / NumPoints;

monte(N,InCircle,NumPoints)->
    X = rand:uniform(),
    Y = rand:uniform(),
    monte(N-1, within(X,Y,InCircle), NumPoints + 1).

within(X,Y,IN)->
          if X*X + Y*Y < 1 -> IN + 1;
          true -> IN
          end.

main(N) -> io:format("PI: ~w~n", [ monte(N) ]).

Output:

Xcoord
6> [monte2:main(X) || X <- [10000000,1000000,100000,10000] ].
PI: 3.1424172
PI: 3.140544
PI: 3.14296
PI: 3.1252
[ok,ok,ok,ok]

ERRE

PROGRAM RANDOM_PI

!
! for rosettacode.org
!

!$DOUBLE

PROCEDURE MONTECARLO(T->RES)
      LOCAL I,N
      FOR I=1 TO T DO
        IF RND(1)^2+RND(1)^2<1 THEN N+=1 END IF
      END FOR
      RES=4*N/T
END PROCEDURE

BEGIN
      RANDOMIZE(TIMER) ! init rnd number generator
      MONTECARLO(1000->RES)     PRINT(RES)
      MONTECARLO(10000->RES)    PRINT(RES)
      MONTECARLO(100000->RES)   PRINT(RES)
      MONTECARLO(1000000->RES)  PRINT(RES)
      MONTECARLO(10000000->RES) PRINT(RES)
END PROGRAM

Output:

 3.136
 3.1468
 3.14392
 3.143824
 3.141514

Euler Math Toolbox

>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):

F#

There is some support and test expressions.

let print x = printfn "%A" x

let MonteCarloPiGreco niter =
    let eng = System.Random()
    let action () =
        let x: float = eng.NextDouble()
        let y: float = eng.NextDouble()
        let res: float = System.Math.Sqrt(x**2.0 + y**2.0)
        if res < 1.0 then
            1
        else
            0
    let res = [ for x in 1..niter do yield action() ]
    let tmp: float = float(List.reduce (+) res) / float(res.Length)
    4.0*tmp

MonteCarloPiGreco 1000 |> print
MonteCarloPiGreco 10000 |> print
MonteCarloPiGreco 100000 |> print

Output:

3.164
3.122
3.1436

Factor

Since Factor lets the user choose the range of the random generator, we use 2^32.

USING: kernel math math.functions random sequences ;

: limit ( -- n ) 2 32 ^ ; inline
: in-circle ( x y -- ? ) limit [ sq ] tri@ [ + ] [ <= ] bi* ;
: rand ( -- r ) limit random ;
: pi ( n -- pi ) [ [ drop rand rand in-circle ] count ] keep / 4 * >float ;

Example use:

10000 pi .
3.1412

Fantom

class MontyCarlo
{
  // assume square/circle of width 1 unit
  static Float findPi (Int samples)
  {
    Int insideCircle := 0
    samples.times 
    {
      x := Float.random
      y := Float.random
      if ((x*x + y*y).sqrt <= 1.0f) insideCircle += 1
    }
    return insideCircle * 4.0f / samples
  }

  public static Void main () 
  {
    [100, 1000, 10000, 1000000, 10000000].each |sample|
    {
      echo ("Sample size $sample gives PI as ${findPi(sample)}")
    }
  }
}

Output:

Sample size 100 gives PI as 3.2
Sample size 1000 gives PI as 3.132
Sample size 10000 gives PI as 3.1612
Sample size 1000000 gives PI as 3.139316
Sample size 10000000 gives PI as 3.1409272

Forth

Works with: GNU Forth

include random.fs

10000 value r

: hit? ( -- ? )
  r random dup *
  r random dup * +
  r dup * < ;

: sims ( n -- hits )
  0 swap 0 do hit? if 1+ then loop ;

1000 sims 4 * . 3232  ok
10000 sims 4 * . 31448  ok
100000 sims 4 * . 313704  ok
1000000 sims 4 * . 3141224  ok
10000000 sims 4 * . 31409400  ok

Fortran

Works with: Fortran version 90 and later

MODULE Simulation
 
   IMPLICIT NONE
 
   CONTAINS
 
   FUNCTION Pi(samples)
     REAL :: Pi
     REAL :: coords(2), length
     INTEGER :: i, in_circle, samples
  
     in_circle = 0
     DO i=1, samples
       CALL RANDOM_NUMBER(coords)
       coords = coords * 2 - 1
       length = SQRT(coords(1)*coords(1) + coords(2)*coords(2))
       IF (length <= 1) in_circle = in_circle + 1
     END DO
     Pi = 4.0 * REAL(in_circle) / REAL(samples)
   END FUNCTION Pi
 
 END MODULE Simulation
  
 PROGRAM MONTE_CARLO
 
   USE Simulation 
   
   INTEGER :: n = 10000
 
   DO WHILE (n <= 100000000)
     WRITE (*,*) n, Pi(n)
     n = n * 10
   END DO
     
 END PROGRAM MONTE_CARLO

Output:

        10000     3.12120
       100000     3.13772
      1000000     3.13934
     10000000     3.14114
    100000000     3.14147

Works with: Fortran version 2008 and later

        program mc
        integer :: n,i
        real(8) :: pi
        n=10000
        do i=1,5
          print*,n,pi(n)
          n = n * 10
        end do
        end program

        function  pi(n)
        integer :: n
        real(8) :: x(2,n),pi
        call random_number(x)
        pi = 4.d0 * dble( count( hypot(x(1,:),x(2,:)) <= 1.d0 ) ) / n
        end function

Futhark

Since Futhark is a pure language, random numbers are implemented using Sobol sequences.

import "futlib/math"

default(f32)

fun dirvcts(): [2][30]i32 =
    [
            [
                536870912, 268435456, 134217728, 67108864, 33554432, 16777216, 8388608, 4194304, 2097152, 1048576, 524288, 262144, 131072, 65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1
            ],
            [
                536870912, 805306368, 671088640, 1006632960, 570425344, 855638016, 713031680, 1069547520, 538968064, 808452096, 673710080, 1010565120, 572653568, 858980352, 715816960, 1073725440, 536879104, 805318656, 671098880, 1006648320, 570434048, 855651072, 713042560, 1069563840, 538976288, 808464432, 673720360, 1010580540, 572662306, 858993459
            ]
    ]


fun grayCode(x: i32): i32 = (x >> 1) ^ x

----------------------------------------
--- Sobol Generator
----------------------------------------
fun testBit(n: i32, ind: i32): bool =
    let t = (1 << ind) in (n & t) == t

fun xorInds(n: i32) (dir_vs: [num_bits]i32): i32 =
    let reldv_vals = zipWith (\ dv i  ->
                                if testBit(grayCode n,i)
                                then dv else 0)
                             dir_vs (iota num_bits)
    in reduce (^) 0 reldv_vals

fun sobolIndI (dir_vs: [m][num_bits]i32, n: i32): [m]i32 =
    map (xorInds n) dir_vs

fun sobolIndR(dir_vs:  [m][num_bits]i32) (n: i32 ): [m]f32 =
    let divisor = 2.0 ** f32(num_bits)
    let arri    = sobolIndI( dir_vs, n )
    in map (\ (x: i32): f32  -> f32(x) / divisor) arri

fun main(n: i32): f32 =
    let rand_nums = map (sobolIndR (dirvcts())) (iota n)
    let dists     = map (\xy ->
                           let (x,y) = (xy[0],xy[1]) in f32.sqrt(x*x + y*y))
                        rand_nums

    let bs        = map (\d -> if d <= 1.0f32 then 1 else 0) dists

    let inside    = reduce (+) 0 bs
    in 4.0f32*f32(inside)/f32(n)

Go

Using standard library math/rand:

package main

import (
    "fmt"
    "math"
    "math/rand"
    "time"
)

func getPi(numThrows int) float64 {
    inCircle := 0
    for i := 0; i < numThrows; i++ {
        //a square with a side of length 2 centered at 0 has 
        //x and y range of -1 to 1
        randX := rand.Float64()*2 - 1 //range -1 to 1
        randY := rand.Float64()*2 - 1 //range -1 to 1
        //distance from (0,0) = sqrt((x-0)^2+(y-0)^2)
        dist := math.Hypot(randX, randY)
        if dist < 1 { //circle with diameter of 2 has radius of 1
            inCircle++
        }
    }
    return 4 * float64(inCircle) / float64(numThrows)
}

func main() {
    rand.Seed(time.Now().UnixNano())
    fmt.Println(getPi(10000))
    fmt.Println(getPi(100000))
    fmt.Println(getPi(1000000))
    fmt.Println(getPi(10000000))
    fmt.Println(getPi(100000000))
}

Output:

3.1164
3.1462
3.142892
3.1419692
3.14149596

Using x/exp/rand:

For very careful Monte Carlo studies, you might consider the subrepository rand library. The random number generator there has some advantages such as better known statistical properties and better use of memory.

package main

import (
    "fmt"
    "math"
    "time"

    "golang.org/x/exp/rand"
)

func getPi(numThrows int) float64 {
    inCircle := 0
    for i := 0; i < numThrows; i++ {
        //a square with a side of length 2 centered at 0 has
        //x and y range of -1 to 1
        randX := rand.Float64()*2 - 1 //range -1 to 1
        randY := rand.Float64()*2 - 1 //range -1 to 1
        //distance from (0,0) = sqrt((x-0)^2+(y-0)^2)
        dist := math.Hypot(randX, randY)
        if dist < 1 { //circle with diameter of 2 has radius of 1
            inCircle++
        }
    }
    return 4 * float64(inCircle) / float64(numThrows)
}

func main() {
    rand.Seed(uint64(time.Now().UnixNano()))
    fmt.Println(getPi(10000))
    fmt.Println(getPi(100000))
    fmt.Println(getPi(1000000))
    fmt.Println(getPi(10000000))
    fmt.Println(getPi(100000000))
}

Haskell

import Control.Monad
import System.Random

getPi throws = do
  results <- replicateM throws one_trial
  return (4 * fromIntegral (sum results) / fromIntegral throws)
  where
    one_trial = do
      rand_x <- randomRIO (-1, 1)
      rand_y <- randomRIO (-1, 1)
      let dist :: Double
          dist = sqrt (rand_x * rand_x + rand_y * rand_y)
      return (if dist < 1 then 1 else 0)

Output:

Example:
 Prelude System.Random Control.Monad> get_pi 10000
 3.1352
 Prelude System.Random Control.Monad> get_pi 100000
 3.15184
 Prelude System.Random Control.Monad> get_pi 1000000
 3.145024

Or, using foldM, and dropping sqrt:

import Control.Monad (foldM, (>=>))
import System.Random (randomRIO)
import Data.Functor ((<&>))

------- APPROXIMATION TO PI BY A MONTE CARLO METHOD ------

monteCarloPi :: Int -> IO Double
monteCarloPi n =
  (/ fromIntegral n) . (4 *) . fromIntegral
    <$> foldM go 0 [1 .. n]
  where
    rnd = randomRIO (0, 1) :: IO Double
    go a _ = rnd >>= ((<&>) rnd . f a)
    f a x y
      | 1 > x ** 2 + y ** 2 = succ a
      | otherwise = a

--------------------------- TEST -------------------------
main :: IO ()
main =
  mapM_
    (monteCarloPi >=> print)
    [1000, 10000, 100000, 1000000]

Output:

For example:

3.244
3.1116
3.14116
3.141396

HicEst

FUNCTION Pi(samples)
   inside = 0
   DO i = 1, samples
      inside = inside + ( (RAN(1)^2 + RAN(1)^2)^0.5 <= 1)
   ENDDO
   Pi = 4 * inside / samples
END

   WRITE(ClipBoard) Pi(1E4) ! 3.1504
   WRITE(ClipBoard) Pi(1E5) ! 3.14204
   WRITE(ClipBoard) Pi(1E6) ! 3.141672
   WRITE(ClipBoard) Pi(1E7) ! 3.1412856

Icon and Unicon

procedure main()
  every t := 10 ^ ( 5 to 9 ) do
     printf("Rounds=%d Pi ~ %r\n",t,getPi(t))
end

link printf

procedure getPi(rounds)
   incircle := 0. 
   every 1 to rounds do 
      if 1 > sqrt((?0 * 2 - 1) ^ 2 + (?0 * 2 - 1) ^ 2) then 
         incircle +:= 1
   return 4 * incircle / rounds
end

Library: Icon Programming Library

printf.icn provides printf

Output:

Rounds=100000 Pi ~ 3.143400
Rounds=1000000 Pi ~ 3.141656
Rounds=10000000 Pi ~ 3.140437
Rounds=100000000 Pi ~ 3.141375
Rounds=1000000000 Pi ~ 3.141604

J

Explicit Solution:

piMC=: monad define "0
  4* y%~ +/ 1>: %: +/ *: <: +: (2,y) ?@$ 0
)

Tacit Solution:

piMCt=: (0.25&* %~ +/@(1 >: [: +/&.:*: _1 2 p. 0 ?@$~ 2&,))"0

Examples:

   piMC 1e6
3.1426
   piMC 10^i.7
4 2.8 3.24 3.168 3.1432 3.14256 3.14014

Alternative Tacit Solution:

pimct=. (4 * +/ % #)@:(1 >: |)@:(? j. ?)@:($&0)"0
   
   (,. pimct) 10 ^ 3 + i.6
  1000   3.168
 10000   3.122
100000 3.13596
   1e6  3.1428
   1e7 3.14158
   1e8 3.14154

Java

public class MC {
	public static void main(String[] args) {
		System.out.println(getPi(10000));
		System.out.println(getPi(100000));
		System.out.println(getPi(1000000));
		System.out.println(getPi(10000000));
		System.out.println(getPi(100000000));
		
	}
	public static double getPi(int numThrows){
		int inCircle= 0;
		for(int i= 0;i < numThrows;i++){
			//a square with a side of length 2 centered at 0 has 
			//x and y range of -1 to 1
			double randX= (Math.random() * 2) - 1;//range -1 to 1
			double randY= (Math.random() * 2) - 1;//range -1 to 1
			//distance from (0,0) = sqrt((x-0)^2+(y-0)^2)
			double dist= Math.sqrt(randX * randX + randY * randY);
			//^ or in Java 1.5+: double dist= Math.hypot(randX, randY);
			if(dist < 1){//circle with diameter of 2 has radius of 1
				inCircle++;
			}
		}
		return 4.0 * inCircle / numThrows;
	}
}

Output:

3.1396
3.14256
3.141516
3.1418692
3.14168604

Works with: Java version 8+

package montecarlo;

import java.util.stream.IntStream;
import java.util.stream.DoubleStream;

import static java.lang.Math.random;
import static java.lang.Math.hypot;
import static java.lang.System.out;

public interface MonteCarlo {
  public static void main(String... arguments) {
    IntStream.of(
      10000,
      100000,
      1000000,
      10000000,
      100000000
    )
      .mapToDouble(MonteCarlo::pi)
      .forEach(out::println)
    ;
  }

  public static double range() {
    //a square with a side of length 2 centered at 0 has 
    //x and y range of -1 to 1
    return (random() * 2) - 1;
  }

  public static double pi(int numThrows){
    long inCircle = DoubleStream.generate(
      //distance from (0,0) = hypot(x, y)
      () -> hypot(range(), range())
    )
      .limit(numThrows)
      .unordered()
      .parallel()
      //circle with diameter of 2 has radius of 1
      .filter(d -> d < 1)
      .count()
    ;
    return (4.0 * inCircle) / numThrows;
  }
}

Output:

3.1556
3.14416
3.14098
3.1419512
3.14160312

JavaScript

ES5

function mcpi(n) {
    var x, y, m = 0;

    for (var i = 0; i < n; i += 1) {
        x = Math.random();
        y = Math.random();

        if (x * x + y * y < 1) {
            m += 1;
        }
    }

    return 4 * m / n;
}

console.log(mcpi(1000));
console.log(mcpi(10000));
console.log(mcpi(100000));
console.log(mcpi(1000000));
console.log(mcpi(10000000));

3.168
3.1396
3.13692
3.140512
3.1417656

ES6

(() => {
    "use strict";

    // --- APPROXIMATION OF PI BY A MONTE CARLO METHOD ---

    // monteCarloPi :: Int -> Float
    const monteCarloPi = n =>
        4 * enumFromTo(1)(n).reduce(a => {
            const [x, y] = [rnd(), rnd()];

            return (x ** 2) + (y ** 2) < 1 ? (
                1 + a
            ) : a;
        }, 0) / n;


    // --------------------- GENERIC ---------------------

    // enumFromTo :: Int -> Int -> [Int]
    const enumFromTo = m =>
        n => Array.from({
            length: 1 + n - m
        }, (_, i) => m + i);


    // rnd :: () -> Float
    const rnd = Math.random;


    // ---------------------- TEST -----------------------
    // From 1000 samples to 10E7 samples
    return enumFromTo(3)(7).forEach(x => {
        const nSamples = 10 ** x;

        console.log(
            `${nSamples} samples: ${monteCarloPi(nSamples)}`
        );
    });
})();

Output:

For example

1000 samples: 3.064
10000 samples: 3.1416
100000 samples: 3.14756
1000000 samples: 3.142536
10000000 samples: 3.142808

jq

Adapted from Wren

Works with: jq

Works with gojq, the Go implementation of jq

jq does not have a built-in PRNG so we will use /dev/urandom as a source of entropy by invoking jq as follows:

# In case gojq is used, trim leading 0s:
function prng {
    cat /dev/urandom | tr -cd '0-9' | fold -w 10 | sed 's/^0*\(.*\)*\(.\)*$/\1\2/'
}    

prng | jq -nMr -f program.jq

program.jq

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

def percent: "\(100000 * . | round / 1000)%";

def pi: 4* (1|atan);

def rfloat: input/1E10;
 
def mcPi:
  . as $n
  | reduce range(0; $n) as $i (0;
      rfloat as $x
      | rfloat as $y
      | if ($x*$x + $y*$y <= 1) then . + 1 else . end)
   | 4 * . / $n ;

"Iterations -> Approx Pi  -> Error",
"----------    ----------    ------",
( pi as $pi
 | range(1; 7)
 | pow(10;.) as $p
 | ($p | mcPi) as $mcpi
 | ((($pi - $mcpi)|length) / $pi) as $error
 | "\($p|lpad(10))    \($mcpi|lpad(10))    \($error|percent|lpad(6))" )

Output:

Iterations -> Approx Pi  -> Error
----------    ----------    ------
        10           2.8    10.873%
       100          3.28    4.406%
      1000         3.172    0.968%
     10000        3.1456    0.128%
    100000       3.13316    0.268%
   1000000      3.139956    0.052%

Jsish

From Javascript ES5 entry, with PRNG seeded during unit testing for reproducibility.

/* Monte Carlo methods, in Jsish */
function mcpi(n) {
    var x, y, m = 0;

    for (var i = 0; i < n; i += 1) {
        x = Math.random();
        y = Math.random();

        if (x * x + y * y < 1) {
            m += 1;
        }
    }

    return 4 * m / n;
}

if (Interp.conf('unitTest')) {
    Math.srand(0);
;    mcpi(1000);
;    mcpi(10000);
;    mcpi(100000);
;    mcpi(1000000);
}

/*
=!EXPECTSTART!=
mcpi(1000) ==> 3.108
mcpi(10000) ==> 3.1236
mcpi(100000) ==> 3.13732
mcpi(1000000) ==> 3.142124
=!EXPECTEND!=
*/

Output:

prompt$ jsish -u monteCarlos.jsi
[PASS] monteCarlos.jsi

Julia

using Printf

function monteπ(n)
    s = count(rand() ^ 2 + rand() ^ 2 < 1 for _ in 1:n)
    return 4s / n
end

for n in 10 .^ (3:8)
    p = monteπ(n)
    println("$(lpad(n, 9)): π ≈ $(lpad(p, 10)), pct.err = ", @sprintf("%2.5f%%", 100 * abs(p - π) / π))
end

Output:

     1000: π ≈      3.224, pct.err = 0.02623%
    10000: π ≈     3.1336, pct.err = 0.254%
   100000: π ≈    3.13468, pct.err = 0.220%
  1000000: π ≈    3.14156, pct.err = 0.001%
 10000000: π ≈  3.1412348, pct.err = 0.011%
100000000: π ≈ 3.14123216, pct.err = 0.011%

K

   sim:{4*(+/{~1<+/(2_draw 0)^2}'!x)%x}

   sim 10000
3.103

   sim'10^!8
4 2.8 3.4 3.072 3.1212 3.14104 3.14366 3.1413

Kotlin

// version 1.1.0

fun mcPi(n: Int): Double {
    var inside = 0
    (1..n).forEach {
        val x = Math.random()
        val y = Math.random()
        if (x * x + y * y <= 1.0) inside++
    }
    return 4.0 * inside / n
}

fun main(args: Array<String>) {   
    println("Iterations -> Approx Pi  -> Error%")
    println("----------    ----------    ------")
    var n = 1_000
    while (n <= 100_000_000) {
        val pi = mcPi(n)
        val err = Math.abs(Math.PI - pi) / Math.PI * 100.0
        println(String.format("%9d  -> %10.8f -> %6.4f", n, pi, err))
        n *= 10
    }
}

Sample output:

Output:

Iterations -> Approx Pi  -> Error%
----------    ----------    ------
     1000  -> 3.12800000 -> 0.4327
    10000  -> 3.15040000 -> 0.2803
   100000  -> 3.14468000 -> 0.0983
  1000000  -> 3.13982000 -> 0.0564
 10000000  -> 3.14182040 -> 0.0072
100000000  -> 3.14160244 -> 0.0003    

Logo

to square :n
  output :n * :n
end
to trial :r
  output less? sum square random :r square random :r  square :r
end
to sim :n :r
  make "hits 0
  repeat :n [if trial :r [make "hits :hits + 1]]
  output 4 * :hits / :n
end

show sim    1000 10000  ; 3.18
show sim   10000 10000  ; 3.1612
show sim  100000 10000  ; 3.145
show sim 1000000 10000  ; 3.140828

LSL

To test it yourself; rez a box on the ground, and add the following as a New Script. (Be prepared to wait... LSL can be slow, but the Servers are typically running thousands of scripts in parallel so what do you expect?)

integer iMIN_SAMPLE_POWER = 0;
integer iMAX_SAMPLE_POWER = 6;
default {
	state_entry() {
		llOwnerSay("Estimating Pi ("+(string)PI+")");
		integer iSample = 0;
		for(iSample=iMIN_SAMPLE_POWER ; iSample<=iMAX_SAMPLE_POWER  ; iSample++) {
			integer iInCircle = 0;
			integer x = 0;
			integer iMaxSamples = (integer)llPow(10, iSample);
			for(x=0 ; x<iMaxSamples ; x++) {
				if(llSqrt(llPow(llFrand(2.0)-1.0, 2.0)+llPow(llFrand(2.0)-1.0, 2.0))<1.0) {
					iInCircle++;
				}
			}
			float fPi = ((4.0*iInCircle)/llPow(10, iSample));
			float fError = llFabs(100.0*(PI-fPi)/PI);
			llOwnerSay((string)iSample+": "+(string)iMaxSamples+" = "+(string)fPi+", Error = "+(string)fError+"%");
		}
		llOwnerSay("Done.");
	}
}

Output:

Estimating Pi (3.141593)
0: 1 = 4.000000, Error = 27.323954%
1: 10 = 4.000000, Error = 27.323954%
2: 100 = 2.880000, Error = 8.326753%
3: 1000 = 3.188000, Error = 1.477192%
4: 10000 = 3.133600, Error = 0.254414%
5: 100000 = 3.138840, Error = 0.087620%
6: 1000000 = 3.142684, Error = 0.034739%
Done.

Lua

function MonteCarlo ( n_throws )
    math.randomseed( os.time() )

    n_inside = 0
    for i = 1, n_throws do
    	if math.random()^2 + math.random()^2 <= 1.0 then
            n_inside = n_inside + 1
    	end
    end    

    return 4 * n_inside / n_throws
end

print( MonteCarlo( 10000 ) )
print( MonteCarlo( 100000 ) )
print( MonteCarlo( 1000000 ) )
print( MonteCarlo( 10000000 ) )

Output:

3.1436
3.13636
3.14376
3.1420188

Mathematica/Wolfram Language

We define a function with variable sample size:

MonteCarloPi[samplesize_Integer] := N[4Mean[If[# > 1, 0, 1] & /@ Norm /@ RandomReal[1, {samplesize, 2}]]]

Example (samplesize=10,100,1000,....10000000):

{#, MonteCarloPi[#]} & /@ (10^Range[1, 7]) // Grid

gives back:

10		3.2
100		3.16
1000		3.152
10000		3.1228
100000		3.14872
1000000		3.1408
10000000	3.14134

monteCarloPi = 4. Mean[UnitStep[1 - Total[RandomReal[1, {2, #}]^2]]] &;
monteCarloPi /@ (10^Range@6)

A less elegant way to solve the problem, is to imagine a (well-trained) monkey, throwing a number of darts at a dartboard.

The darts land randomly on the board, at different x and y coordinates. We want to know how many darts land inside the circle. We then guess Pi by dividing the total number of darts inside the circle by the total number of darts thrown (assuming they all hit the square board), and multiplying the whole lot by 4.

We create a function MonkeyDartsPi, which can take a variable number of throws as input:

MonkeyDartsPi[numberOfThrows_] := (
xyCoordinates = RandomReal[{0, 1}, {numberOfThrows, 2}];
InsideCircle = Length[Select[Total[xyCoordinates^2, {2}],#<=1&]] ;
4*N[InsideCircle / Length[xyCoordinates],1+Log10[numberOfThrows]])

We do several runs with a larger number of throws each time, increasing by powers of 10.

Grid[Table[{n, MonkeyDartsPi[n]}, {n, 10^Range[7]} ], Alignment -> Left]

We see that as the number of throws increases, we get closer to the value of Pi:

10 		2.8
100 		3.20
1000 		3.176
10000 		3.1356
100000 		3.13700
1000000 	3.142624
10000000 	3.1416328

MATLAB

See: Monte Carlo Simulation in MATLAB for more examples

The first example given is not vectorized. MATLAB has a self-imposed memory limitation that prevents this simulation from having more than 3 decimal digits of accuracy. Because of this limitation it is best to vectorize the code as much as possible so extra memory isn't consumed by unneeded variables. Therefore, I have provided a second solution that is maximally vectorized.

Minimally Vectorized:

function piEstimate = monteCarloPi(numDarts)

    %The square has a sides of length 2, which means the circle has radius
    %1.
    
    %Generate a table of random x-y value pairs in the range [0,1] sampled
    %from the uniform distribution for each axis.
    darts = rand(numDarts,2);
    
    %Any darts that are in the circle will have position vector whose
    %length is less than or equal to 1 squared.
    dartsInside = ( sum(darts.^2,2) <= 1 );
    
    piEstimate = 4*sum(dartsInside)/numDarts;

end

Completely Vectorized:

function piEstimate = monteCarloPi(numDarts)
    
    piEstimate = 4*sum( sum(rand(numDarts,2).^2,2) <= 1 )/numDarts;

end

Output:

>> monteCarloPi(7000000)

ans =

   3.141512000000000

Maxima

load("distrib");
approx_pi(n):= block(
  [x: random_continuous_uniform(0, 1, n),
   y: random_continuous_uniform(0, 1, n),
   r, cin: 0, listarith: true],
   r: x^2 + y^2,
   for r0 in r do if r0<1 then cin: cin + 1,
   4*cin/n);
 
float(approx_pi(100));

MAXScript

fn monteCarlo iterations =
(
    radius = 1.0
    pointsInCircle = 0
    for i in 1 to iterations do
    (
        testPoint = [(random -radius radius), (random -radius radius)]
        if length testPoint <= radius then
        (
            pointsInCircle += 1
        )
    )
    4.0 * pointsInCircle / iterations
)

МК-61/52

П0	П1	0	П4	СЧ	x^2	^	СЧ	x^2	+
1	-	x<0	15	КИП4	L0	04	ИП4	4	*
ИП1	/	С/П

Example: for n = 1000 the output is 3.152.

Nim

import math, random

randomize()
 
proc pi(nthrows: float): float =
  var inside = 0.0
  for i in 1..int64(nthrows):
    if hypot(rand(1.0), rand(1.0)) < 1:
      inside += 1
  result = 4 * inside / nthrows
 
for n in [10e4, 10e6, 10e7, 10e8]:
  echo pi(n)

Output:

3.15336
3.1405116
3.14163332
3.141486144

OCaml

let get_pi throws =
  let rec helper i count =
    if i = throws then count
    else
      let rand_x = Random.float 2.0 -. 1.0
      and rand_y = Random.float 2.0 -. 1.0 in
      let dist = sqrt (rand_x *. rand_x +. rand_y *. rand_y) in
      if dist < 1.0 then
        helper (i+1) (count+1)
      else
        helper (i+1) count
  in float (4 * helper 0 0) /. float throws

Example:

# get_pi 10000;;
- : float = 3.15
# get_pi 100000;;
- : float = 3.13272
# get_pi 1000000;;
- : float = 3.143808
# get_pi 10000000;;
- : float = 3.1421704
# get_pi 100000000;;
- : float = 3.14153872

Octave

function p = montepi(samples)
  in_circle = 0;
  for samp = 1:samples
    v = [ unifrnd(-1,1), unifrnd(-1,1) ];
    if ( v*v.' <= 1.0 )
      in_circle++;
    endif
  endfor
  p = 4*in_circle/samples;
endfunction

l = 1e4;
while (l < 1e7)
  disp(montepi(l));
  l *= 10;
endwhile

Since it runs slow, I've stopped it at the second iteration, obtaining:

 3.1560
 3.1496

Much faster implementation

function result = montepi(n)
  result = sum(rand(1,n).^2+rand(1,n).^2<1)/n*4;
endfunction

PARI/GP

MonteCarloPi(tests)=4.*sum(i=1,tests,norml2([random(1.),random(1.)])<1)/tests;

A hundred million tests (about a minute) yielded 3.14149000, slightly more precise (and round!) than would have been expected. A million gave 3.14162000 and a thousand 3.14800000.

Pascal

Library: Math

Program MonteCarlo(output);

uses
  Math;

function MC_Pi(expo: integer): real;
  var
    x, y: real;
    i, hits, samples: longint;
  begin
    samples := 10**expo;
    hits := 0;
    randomize;
    for i := 1 to samples do
    begin
      x := random;
      y := random;
      if sqrt(x*x + y*y) < 1.0 then
        inc(hits);
    end;
    MC_Pi := 4.0 * hits / samples;
  end;

var
  i: integer;
begin
  for i := 4 to 8 do
    writeln (10**i, ' samples give ', MC_Pi(i):7:5, ' as pi.');
end.

Output:

:> ./MonteCarlo
10000 samples give 3.14480 as pi.
100000 samples give 3.14484 as pi.
1000000 samples give 3.13970 as pi.
10000000 samples give 3.14100 as pi.
100000000 samples give 3.14162 as pi.

Perl

sub pi {
  my $nthrows = shift;
  my $inside = 0;
  foreach (1 .. $nthrows) {
    my $x = rand() * 2 - 1;
    my $y = rand() * 2 - 1;
    if (sqrt($x*$x + $y*$y) < 1) {
      $inside++;
    }
  }
  return 4 * $inside / $nthrows;
}

printf "%9d: %07f\n", $_, pi($_) for 10**4, 10**6;

Output:

    10000: 3.132000
  1000000: 3.141596

Phix

with javascript_semantics
integer N = 100
for i=1 to 6 do
    integer inside = 0
    for n=1 to N do
        integer x = rand(N),
                y = rand(N)
        inside += (x*x+y*y<N*N)
    end for
    pp({N,4*inside/N})
    N *= 10
end for

Output:

{100,3.2}
{1000,3.116}
{10000,3.1736}
{100000,3.13996}
{1000000,3.141856}
{10000000,3.1415728}

PHP

<?
$loop = 1000000; # loop to 1,000,000
$count = 0;
for ($i=0; $i<$loop; $i++) {
  $x = rand() / getrandmax();
  $y = rand() / getrandmax();
  if(($x*$x) + ($y*$y)<=1) $count++;
}
echo "loop=".number_format($loop).", count=".number_format($count).", pi=".($count/$loop*4);
?>

Output:

loop=1,000,000, count=785,462, pi=3.141848

Picat

Some general Monte Carlo simulators. N is the number of runs, F is the simulation function.

Using while loop

sim1(N, F) = C =>
  C = 0,
  I = 0,
  while (I <= N) 
    C := C + apply(F),
    I := I + 1
  end.

List comprehension

This is simpler, but slightly slower than using while loop.

sim2(N, F) = sum([apply(F) : _I in 1..N]).

Recursion

sim_rec(N,F) = S => 
  sim_rec(N,N,F,0,S).
sim_rec(0,_N,_F,S,S).
sim_rec(C,N,F,S0,S) :-
  S1 = S0 + apply(F),
  sim_rec(C-1,N,F,S1,S).

Test

Of the three different MC simulators, sim_rec/2 (using recursion) is slightly faster than the other two (sim1/2 and sim2/2) which have about the same speed.

go =>
   foreach(N in 0..7)
     sim_pi(10**N)
   end,
   nl.

% The specific pi simulation
sim_pi(N) =>
  Inside = sim(N,pi_f),
  MyPi = 4.0*Inside/N,
  Pi = math.pi,
  println([n=N, myPi=MyPi, diff=Pi-MyPi]).

% The simulation function: 
%   returns 1 if success, 0 otherwise
pi_f() = cond(frand()**2 + frand()**2 <= 1, 1, 0).

Output:

[n = 1,myPi = 4.0,diff = -0.858407346410207]
[n = 10,myPi = 3.2,diff = -0.058407346410207]
[n = 100,myPi = 3.12,diff = 0.021592653589793]
[n = 1000,myPi = 3.152,diff = -0.010407346410207]
[n = 10000,myPi = 3.1672,diff = -0.025607346410207]
[n = 100000,myPi = 3.13888,diff = 0.002712653589793]
[n = 1000000,myPi = 3.14192,diff = -0.000327346410207]
[n = 10000000,myPi = 3.1408988,diff = 0.000693853589793]

PicoLisp

(de carloPi (Scl)
   (let (Dim (** 10 Scl)  Dim2 (* Dim Dim)  Pi 0)
      (do (* 4 Dim)
         (let (X (rand 0 Dim)  Y (rand 0 Dim))
            (when (>= Dim2 (+ (* X X) (* Y Y)))
               (inc 'Pi) ) ) )
      (format Pi Scl) ) )

(for N 6
   (prinl (carloPi N)) )

Output:

3.4
3.23
3.137
3.1299
3.14360
3.140964

PowerShell

Works with: PowerShell version 2

function Get-Pi ($Iterations = 10000) {
    $InCircle = 0
    for ($i = 0; $i -lt $Iterations; $i++) {
        $x = Get-Random 1.0
        $y = Get-Random 1.0
        if ([Math]::Sqrt($x * $x + $y * $y) -le 1) {
            $InCircle++
        }
    }
    $Pi = [decimal] $InCircle / $Iterations * 4
    $RealPi = [decimal] "3.141592653589793238462643383280"
    $Diff = [Math]::Abs(($Pi - $RealPi) / $RealPi * 100)
    New-Object PSObject `
        | Add-Member -PassThru NoteProperty Iterations $Iterations `
        | Add-Member -PassThru NoteProperty Pi $Pi `
        | Add-Member -PassThru NoteProperty "% Difference" $Diff
}

This returns a custom object with appropriate properties which automatically enables a nice tabular display:

PS Home:\> 10,100,1e3,1e4,1e5,1e6 | ForEach-Object { Get-Pi $_ }

 Iterations          Pi                      % Difference
 ----------          --                      ------------
         10         3,6    14,591559026164641753596309630
        100        3,40     8,225361302488828322840959090
       1000       3,208    2,1138114877600474293158225800
      10000      3,1444    0,0893606116311387583356211100
     100000     3,14712    0,1759409006731298209938938800
    1000000    3,141364    0,0072782698142600895432451100

Python

At the interactive prompt

Python 2.6rc2 (r26rc2:66507, Sep 18 2008, 14:27:33) [MSC v.1500 32 bit (Intel)] on win32 IDLE 2.6rc2

One use of the "sum" function is to count how many times something is true (because True = 1, False = 0):

>>> import random, math
>>> throws = 1000
>>> 4.0 * sum(math.hypot(*[random.random()*2-1
	                 for q in [0,1]]) < 1
              for p in xrange(throws)) / float(throws)
3.1520000000000001
>>> throws = 1000000
>>> 4.0 * sum(math.hypot(*[random.random()*2-1
	                 for q in [0,1]]) < 1
              for p in xrange(throws)) / float(throws)
3.1396359999999999
>>> throws = 100000000
>>> 4.0 * sum(math.hypot(*[random.random()*2-1
	                 for q in [0,1]]) < 1
              for p in xrange(throws)) / float(throws)
3.1415666400000002

As a program using a function

from random import random
from math import hypot
try:
    import psyco
    psyco.full()
except:
    pass

def pi(nthrows):
    inside = 0
    for i in xrange(nthrows):
        if hypot(random(), random()) < 1:
            inside += 1
    return 4.0 * inside / nthrows

for n in [10**4, 10**6, 10**7, 10**8]:
    print "%9d: %07f" % (n, pi(n))

Faster implementation using Numpy

import numpy as np

n = input('Number of samples: ')
print np.sum(np.random.rand(n)**2+np.random.rand(n)**2<1)/float(n)*4

Quackery

Translation of: Forth

  [ $ "bigrat.qky" loadfile ] now!

  [ [ 64 bit ] constant
    dup random dup *
    over random dup * +
    swap dup * < ]            is hit    (   --> b )

  [ 0 swap times
      [ hit if 1+ ] ]         is sims   ( n --> n )

  [ dup echo say " trials "
    dup sims 4 *
    swap 20 point$ echo$ cr ] is trials ( n -->   )

' [ 10 100 1000 10000 100000 1000000 ] witheach trials

Output:

10 trials 2.8
100 trials 3.2
1000 trials 3.172
10000 trials 3.1484
100000 trials 3.1476
1000000 trials 3.142256

R

# nice but not suitable for big samples!
monteCarloPi <- function(samples) {
  x <- runif(samples, -1, 1) # for big samples, you need a lot of memory!
  y <- runif(samples, -1, 1)
  l <- sqrt(x*x + y*y)
  return(4*sum(l<=1)/samples)
}

# this second function changes the samples number to be
# multiple of group parameter (default 100).
monteCarlo2Pi <- function(samples, group=100) {
  lim <- ceiling(samples/group)
  olim <- lim
  c <- 0
  while(lim > 0) {
    x <- runif(group, -1, 1)
    y <- runif(group, -1, 1)
    l <- sqrt(x*x + y*y)
    c <- c + sum(l <= 1)
    lim <- lim - 1
  }
  return(4*c/(olim*group))
}

print(monteCarloPi(1e4))
print(monteCarloPi(1e5))
print(monteCarlo2Pi(1e7))

Racket

#lang racket

(define (in-unit-circle? x y) (<= (sqrt (+ (sqr x) (sqr y))) 1))
;; point in ([-1,1], [-1,1])
(define (random-point-in-2x2-square) (values (* 2 (- (random) 1/2)) (* 2 (- (random) 1/2))))

;; Area of circle is (pi r^2). r is 1, area of circle is pi
;; Area of square is 2^2 = 4
;; There is a pi/4 chance of landing in circle
;; .: pi = 4*(proportion passed) = 4*(passed/samples)
(define (passed:samples->pi passed samples) (* 4 (/ passed samples)))

;; generic kind of monte-carlo simulation
(define (monte-carlo run-length report-frequency
                     sample-generator pass?
                     interpret-result)
  (let inner ((samples 0) (passed 0) (cnt report-frequency))
    (cond
      [(= samples run-length) (interpret-result passed samples)]
      [(zero? cnt) ; intermediate report
       (printf "~a samples of ~a: ~a passed -> ~a~%"
               samples run-length passed (interpret-result passed samples))
       (inner samples passed report-frequency)]
      [else
       (inner (add1 samples)
              (if (call-with-values sample-generator pass?)
                  (add1 passed) passed) (sub1 cnt))])))

;; (monte-carlo ...) gives an "exact" result... which will be a fraction.
;; to see how it looks as a decimal we can exact->inexact it
(let ((mc (monte-carlo 10000000 1000000 random-point-in-2x2-square in-unit-circle? passed:samples->pi)))
  (printf "exact = ~a~%inexact = ~a~%(pi - guess) = ~a~%" mc (exact->inexact mc) (- pi mc)))

Output:

1000000 samples of 10000000: 785763 passed -> 785763/250000
2000000 samples of 10000000: 1571487 passed -> 1571487/500000
3000000 samples of 10000000: 2356776 passed -> 98199/31250
4000000 samples of 10000000: 3141924 passed -> 785481/250000
5000000 samples of 10000000: 3927540 passed -> 196377/62500
6000000 samples of 10000000: 4713072 passed -> 98189/31250
7000000 samples of 10000000: 5498300 passed -> 54983/17500
8000000 samples of 10000000: 6283199 passed -> 6283199/2000000
9000000 samples of 10000000: 7068065 passed -> 1413613/450000
exact = 3926793/1250000
inexact = 3.1414344
(pi - guess) = 0.00015825358979304482

A little more Racket-like is the use of an iterator (in this case for/fold), which is clearer than an inner function:

#lang racket
(define (in-unit-circle? x y) (<= (sqrt (+ (sqr x) (sqr y))) 1))
;; Good idea made in another task that:
;;  The proportions of hits is the same in the unit square and 1/4 of a circle.
;; point in ([0,1], [0,1])
(define (random-point-in-unit-square) (values (random) (random)))
;; generic kind of monte-carlo simulation
;; Area of circle is (pi r^2). r is 1, area of circle is pi
;; Area of square is 2^2 = 4
;; There is a pi/4 chance of landing in circle
;; .: pi = 4*(proportion passed) = 4*(passed/samples)
(define (passed:samples->pi passed samples) (* 4 (/ passed samples)))

(define (monte-carlo/2 run-length report-frequency sample-generator pass? interpret-result)
  (interpret-result
   (for/fold ((pass 0))
     ([n (in-range run-length)]
      #:when (when (and (not (zero? n)) (zero? (modulo n report-frequency)))
               (printf "~a samples of ~a: ~a passed -> ~a~%"
                       n run-length pass (interpret-result pass n)))
      #:when (call-with-values sample-generator pass?))
     (add1 pass))
   run-length))

;; (monte-carlo ...) gives an "exact" result... which will be a fraction.
;; to see how it looks as a decimal we can exact->inexact it
(let ((mc (monte-carlo/2 10000000 1000000 random-point-in-unit-square in-unit-circle? passed:samples->pi)))
  (printf "exact = ~a~%inexact = ~a~%(pi - guess) = ~a~%" mc (exact->inexact mc) (- pi mc)))

[Similar output]

Raku

(formerly Perl 6)

Works with: rakudo version 2015-09-24

We'll consider the upper-right quarter of the unitary disk centered at the origin. Its area is $\pi \over 4$.

my @random_distances = ([+] rand**2 xx 2) xx *;

sub approximate_pi(Int $n) {
    4 * @random_distances[^$n].grep(* < 1) / $n
}

say "Monte-Carlo π approximation:";
say "$_ iterations:  ", approximate_pi $_
    for 100, 1_000, 10_000;

Output:

Monte-Carlo π approximation:
100 iterations:  2.88
1000 iterations:  3.096
10000 iterations:  3.1168

We don't really need to write a function, though. A lazy list would do:

my @pi = ([\+] 4 * (1 > [+] rand**2 xx 2) xx *) Z/ 1 .. *;
say @pi[10, 1000, 10_000];

REXX

A specific─purpose commatizer function is included to format the number of iterations.

/*REXX program computes and displays the value of  pi÷4  using the Monte Carlo algorithm*/
numeric digits 20                                /*use 20 decimal digits to handle args.*/
parse arg times chunk digs r? .                  /*does user want a specific number?    */
if times=='' | times==","  then times=   5e12    /*five trillion should do it, hopefully*/
if chunk=='' | chunk==","  then chunk= 100000    /*perform Monte Carlo in  100k  chunks.*/
if digs =='' |  digs==","  then  digs=     99    /*indicates to use length of  PI - 1.  */
if datatype(r?, 'W')       then call random ,,r? /*Is there a random seed?  Then use it.*/
                                                 /* [↓]  pi meant to line─up with a SAY.*/
           pi= 3.141592653589793238462643383279502884197169399375105820974944592307816406
pi= strip( left(pi, digs + length(.) ) )         /*obtain length of pi to what's wanted.*/
numeric digits length(pi) - 1                    /*define decimal digits as length PI -1*/
say '                    1         2         3         4         5         6         7   '
say 'scale:    1·234567890123456789012345678901234567890123456789012345678901234567890123'
say                                              /* [↑]  a two─line scale for showing pi*/
say 'true pi= '       pi"+"                      /*we might as well brag about true  pi.*/
say                                              /*display a blank line for separation. */
limit   = 10000 - 1                              /*REXX random generates only integers. */
limitSq = limit **2                              /*··· so, instead of one, use limit**2.*/
accuracy= 0                                      /*accuracy of Monte Carlo pi  (so far).*/
@reps= 'repetitions:  Monte Carlo  pi  is'       /*a handy─dandy short literal for a SAY*/
!= 0                                             /*!:  is the accuracy of pi  (so far). */
       do j=1  for times % chunk
                     do chunk                    /*do Monte Carlo, one chunk at─a─time. */
                     if random(, limit)**2 + random(, limit)**2 <= limitSq   then != ! + 1
                     end   /*chunk*/
       reps= chunk * j                           /*calculate the number of repetitions. */
       _= compare(4*! / reps, pi)                /*compare apples and  ···  crabapples. */
       if _<=accuracy  then iterate              /*Not better accuracy?  Keep truckin'. */
       say right(commas(reps), 20) @reps  'accurate to'  _-1  "places."  /*─1≡dec. point*/
       accuracy= _                               /*use this accuracy for next baseline. */
       end   /*j*/
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: procedure; arg _; do k=length(_)-3  to 1  by -3; _=insert(',',_,k); end;  return _

output   when using the default input:

                    1         2         3         4         5         6         7
scale:    1·234567890123456789012345678901234567890123456789012345678901234567890123

true pi=  3.141592653589793238462643383279502884197169399375105820974944592307816406+

             10,000 repetitions:  Monte Carlo  pi  is accurate to 3 places.
             50,000 repetitions:  Monte Carlo  pi  is accurate to 4 places.
            850,000 repetitions:  Monte Carlo  pi  is accurate to 5 places.
            890,000 repetitions:  Monte Carlo  pi  is accurate to 6 places.
          5,130,000 repetitions:  Monte Carlo  pi  is accurate to 7 places.
          8,620,000 repetitions:  Monte Carlo  pi  is accurate to 8 places.
         10,390,000 repetitions:  Monte Carlo  pi  is accurate to 9 places.

For more example runs using REXX,   see the   discussion   page.

Ring

decimals(8)
see "monteCarlo(1000) = " + monteCarlo(1000) + nl
see "monteCarlo(10000) = " + monteCarlo(10000) + nl
see "monteCarlo(100000) = " + monteCarlo(100000) + nl
 
func monteCarlo t
     n=0
     for i = 1 to t
         if sqrt(pow(random(1),2) + pow(random(1),2)) <= 1 n += 1 ok
     next
     t = (4 * n) / t
     return t

Output:

monteCarlo(1000) = 3.11600000
monteCarlo(10000) = 3.00320000
monteCarlo(100000) = 2.99536000

RPL

≪ 0 
   1 3 PICK START
     RAND SQ RAND SQ + 1 < +
   NEXT
   SWAP / 4 * 
≫ 'MCARL' STO

100 MCARL
1000 MCARL
10000 MCARL
100000 MCARL

Output:

4: 3.2
3: 3.084
2: 3.1684
1: 3.14154

Ruby

def approx_pi(throws)
  times_inside = throws.times.count {Math.hypot(rand, rand) <= 1.0}
  4.0 * times_inside / throws
end

[1000, 10_000, 100_000, 1_000_000, 10_000_000].each do |n| 
   puts "%8d samples: PI = %s" % [n, approx_pi(n)]
end

Output:

    1000 samples: PI = 3.2
   10000 samples: PI = 3.14
  100000 samples: PI = 3.13244
 1000000 samples: PI = 3.145124
10000000 samples: PI = 3.1414788

Rust

extern crate rand;

use rand::Rng;
use std::f64::consts::PI;

// `(f32, f32)` would be faster for some RNGs (including `rand::thread_rng` on 32-bit platforms
// and `rand::weak_rng` as of rand v0.4) as `next_u64` combines two `next_u32`s if not natively
// supported by the RNG.  It would less accurate however.
fn is_inside_circle((x, y): (f64, f64)) -> bool {
    x * x + y * y <= 1.0
}

fn simulate<R: Rng>(rng: &mut R, samples: usize) -> f64 {
    let mut count = 0;
    for _ in 0..samples {
        if is_inside_circle(rng.gen()) {
            count += 1;
        }
    }
    (count as f64) / (samples as f64)
}

fn main() {
    let mut rng = rand::weak_rng();

    println!("Real pi: {}", PI);

    for samples in (3..9).map(|e| 10_usize.pow(e)) {
        let estimate = 4.0 * simulate(&mut rng, samples);
        let deviation = 100.0 * (1.0 - estimate / PI).abs();
        println!("{:9}: {:<11} dev: {:.5}%", samples, estimate, deviation);
    }
}

Output:

Real pi: 3.141592653589793
     1000: 3.212       dev: 2.24114%
    10000: 3.156       dev: 0.45860%
   100000: 3.14112     dev: 0.01505%
  1000000: 3.14122     dev: 0.01186%
 10000000: 3.1408112   dev: 0.02487%
100000000: 3.14186092  dev: 0.00854%

Scala

object MonteCarlo {
  private val random = new scala.util.Random

  /** Returns a random number between -1 and 1 */
  def nextThrow: Double = (random.nextDouble * 2.0) - 1.0

  /** Returns true if the argument point would be 'inside' the unit circle with
    * center at the origin, and bounded by a square with side lengths of 2
    * units. */
  def insideCircle(pt: (Double, Double)): Boolean = pt match {
    case (x, y) => (x * x) + (y * y) <= 1.0
  }
  
  /** Runs the simulation the specified number of times. Uses the result to 
    * estimate a value of pi */
  def simulate(times: Int): Double = {
    val inside = Iterator.tabulate (times) (_ => (nextThrow, nextThrow)) count insideCircle
    inside.toDouble / times.toDouble * 4.0
  }

  def main(args: Array[String]): Unit = {
    val sims = Seq(10000, 100000, 1000000, 10000000, 100000000)
    sims.foreach { n =>
      println(n+" simulations; pi estimation: "+ simulate(n))
    }
  }
}

Output:

10000 simulations; pi estimation: 3.1492
100000 simulations; pi estimation: 3.1396
1000000 simulations; pi estimation: 3.14208
10000000 simulations; pi estimation: 3.1409944
100000000 simulations; pi estimation: 3.1414386

Seed7

$ include "seed7_05.s7i";
  include "float.s7i";

const func float: pi (in integer: throws) is func
  result
    var float: pi is 0.0;
  local
    var integer: throw is 0;
    var integer: inside is 0;
  begin
    for throw range 1 to throws do
      if rand(0.0, 1.0) ** 2 + rand(0.0, 1.0) ** 2 <= 1.0 then
        incr(inside);
      end if;
    end for;
    pi := flt(4 * inside) / flt(throws);
  end func;

const proc: main is func
  begin
    writeln("    10000: " <& pi(    10000) digits 5);
    writeln("   100000: " <& pi(   100000) digits 5);
    writeln("  1000000: " <& pi(  1000000) digits 5);
    writeln(" 10000000: " <& pi( 10000000) digits 5);
    writeln("100000000: " <& pi(100000000) digits 5);
  end func;

Output:

    10000: 3.14520
   100000: 3.15000
  1000000: 3.14058
 10000000: 3.14223
100000000: 3.14159

SequenceL

First solution is serial due to the use of random numbers. Will always give the same result for a given n and seed

import <Utilities/Random.sl>;
import <Utilities/Conversion.sl>;

main(args(2)) := monteCarlo(stringToInt(args[1]), stringToInt(args[2]));

monteCarlo(n, seed) :=
	let
		totalHits := monteCarloHelper(n, seedRandom(seed), 0);
	in
		(totalHits / intToFloat(n))*4.0;

monteCarloHelper(n, generator, result) :=
	let
		xRand := getRandom(generator);
		x := xRand.Value/(generator.RandomMax + 1.0);
		yRand := getRandom(xRand.Generator);
		y := yRand.Value/(generator.RandomMax + 1.0);
		
		newResult := result + 1 when x^2 + y^2 < 1.0 else
					 result;
	in
		result when n < 0 else
		monteCarloHelper(n - 1, yRand.Generator, newResult);

The second solution will run in parallel. It will also always give the same result for a given n and seed. (Note, the function monteCarloHelper is the same in both versions).

import <Utilities/Random.sl>;
import <Utilities/Conversion.sl>;

main(args(2)) := monteCarlo(stringToInt(args[1]), stringToInt(args[2]));

chunks := 100;
monteCarlo3(n, seed) :=
	let
		newSeeds := getRandomSequence(seedRandom(seed), chunks).Value;
		totalHits := monteCarloHelper(n / chunks, seedRandom(newSeeds), 0);
	in
		(sum(totalHits) / intToFloat((n / chunks)*chunks))*4.0;

monteCarloHelper(n, generator, result) :=
	let
		xRand := getRandom(generator);
		x := xRand.Value/(generator.RandomMax + 1.0);
		yRand := getRandom(xRand.Generator);
		y := yRand.Value/(generator.RandomMax + 1.0);
		
		newResult := result + 1 when x^2 + y^2 < 1.0 else
					 result;
	in
		result when n < 0 else
		monteCarloHelper(n - 1, yRand.Generator, newResult);

Sidef

func monteCarloPi(nthrows) {
    4 * (^nthrows -> count_by {
        hypot(1.rand(2) - 1, 1.rand(2) - 1) < 1
    }) / nthrows
}

for n in [1e2, 1e3, 1e4, 1e5, 1e6] {
    printf("%9d: %07f\n", n, monteCarloPi(n))
}

Output:

      100: 3.320000
     1000: 3.120000
    10000: 3.169600
   100000: 3.138920
  1000000: 3.142344

SparForte

As a structured script.

#!/usr/local/bin/spar
pragma annotate( summary, "monte" )
              @( description, "A Monte Carlo Simulation is a way of approximating the" )
              @( description, "value of a function where calculating the actual value is" )
              @( description, "difficult or impossible. It uses random sampling to define" )
              @( description, "constraints on the value and then makes a sort of 'best" )
              @( description, "guess.'" )
              @( description, "" )
              @( description, "Write a function to run a simulation like this with a" )
              @( description, "variable number of random points to select. Also, show the" )
              @( description, "results of a few different sample sizes. For software" )
              @( description, "where the number pi is not built-in, we give pi to a couple" )
              @( description, "of digits: 3.141592653589793238462643383280 " )
              @( see_also, "http://rosettacode.org/wiki/Monte_Carlo_methods" )
              @( author, "Ken O. Burtch" );
pragma license( unrestricted );

pragma restriction( no_external_commands );

procedure monte is
 
   function pi_estimate (throws : positive) return float is
      inside : natural := 0;
   begin
      for throw in 1..throws loop
         if numerics.random ** 2 + numerics.random ** 2 <= 1.0 then
            inside := @ + 1;
         end if;
      end loop;
      return 4.0 * float (inside) / float (throws);
   end pi_estimate;

begin
   ? "      1_000:" & strings.image (pi_estimate (      1_000))
   @ "     10_000:" & strings.image (pi_estimate (     10_000))
   @ "    100_000:" & strings.image (pi_estimate (    100_000))
   @ "  1_000_000:" & strings.image (pi_estimate (  1_000_000));
end monte;

Stata

program define mcdisk
	clear all
	quietly set obs `1'
	gen x=2*runiform()
	gen y=2*runiform()
	quietly count if (x-1)^2+(y-1)^2<1
	display 4*r(N)/_N
end

. mcdisk 10000
3.1424

. mcdisk 1000000
3.141904

. mcdisk 100000000
3.1416253

Swift

Translation of: JavaScript

import Foundation

func mcpi(sampleSize size:Int) -> Double {
    var x = 0 as Double
    var y = 0 as Double
    var m = 0 as Double
    
    for i in 0..<size {
        x = Double(arc4random()) / Double(UINT32_MAX)
        y = Double(arc4random()) / Double(UINT32_MAX)
        
        if ((x * x) + (y * y) < 1) {
            m += 1
        }
    }
    
    return (4.0 * m) / Double(size)
}

println(mcpi(sampleSize: 100))
println(mcpi(sampleSize: 1000))
println(mcpi(sampleSize: 10000))
println(mcpi(sampleSize: 100000))
println(mcpi(sampleSize: 1000000))
println(mcpi(sampleSize: 10000000))
println(mcpi(sampleSize: 100000000))

Output:

3.08
3.128
3.1548
3.149
3.142032
3.1414772
3.14166832

Tcl

proc pi {samples} {
    set i 0
    set inside 0
    while {[incr i] <= $samples} {
        if {sqrt(rand()**2 + rand()**2) <= 1.0} {
            incr inside
        }
    }
    return [expr {4.0 * $inside / $samples}]
}

puts "PI is approx [expr {atan(1)*4}]\n"
foreach runs {1e2 1e4 1e6 1e8} {
    puts "$runs => [pi $runs]"
}

result

PI is approx 3.141592653589793

1e2 => 2.92
1e4 => 3.1344
1e6 => 3.141924
1e8 => 3.14167724

Ursala

#import std
#import flo

mcp "n" = times/4. div\float"n" (rep"n" (fleq/.5+ sqrt+ plus+ ~~ sqr+ minus/.5+ rand)?/~& plus/1.) 0.

Here's a walk through.

• mcp "n" = ... defines a function named mcp in terms of a dummy variable "n", which will be the number of iterations used in the simulation
• rand ignores its argument and returns a uniformly distributed number between 0 and 1
• minus/.5 is composed with rand to compute the difference between the random number and 0.5
• sqr squares the difference
• ~~ says to apply the function twice and return the pair of results
• plus composed with that adds the pair of results
• sqrt takes the square root of the sum
• fleq/.5 is floating point comparison with a fixed right side of .5, returning true if its argument is greater or equal
• Everything from fleq to rand forms the predicate for the ? conditional operator.
• If the condition is true, the identity function is applied, ~&
• If the condition is false, the plus/1. function is applied, which adds one to its argument.
• rep"n" applied to a function has the effect of composing that function with itself "n" times, with "n" in this case being the parameter to the mcp function.
• The function being repeated "n" times is applied to an argument of 0.
• A division of the result by the number "n" converted to a floating point value is performed by div\float"n".
• The result of the division is quadrupled by times/4..

test program:

#cast %eL

pis = mcp* <10,100,1000,10000,100000,1000000>

Output:

<
   2.800000e+00,
   3.600000e+00,
   3.164000e+00,
   3.118800e+00,
   3.144480e+00,
   3.141668e+00>

Wren

Translation of: Kotlin

Library: Wren-fmt

import "random" for Random
import "./fmt" for Fmt

var rand = Random.new()

var mcPi = Fn.new { |n|
    var inside = 0
    for (i in 1..n) {
        var x = rand.float()
        var y = rand.float()
        if (x*x + y*y <= 1) inside = inside + 1
    }
    return 4 * inside / n
}

System.print("Iterations -> Approx Pi  -> Error\%")
System.print("----------    ----------    ------")
var n = 1000
while (n <= 1e8) {
    var pi = mcPi.call(n)
    var err = (Num.pi - pi).abs / Num.pi * 100.0
    Fmt.print("$9d  -> $10.8f -> $6.4f", n, pi, err)
    n = n * 10
}

Output:

Sample run:

Iterations -> Approx Pi  -> Error%
----------    ----------    ------
     1000  -> 3.21200000 -> 2.2411
    10000  -> 3.16720000 -> 0.8151
   100000  -> 3.13944000 -> 0.0685
  1000000  -> 3.14048000 -> 0.0354
 10000000  -> 3.14191240 -> 0.0102
100000000  -> 3.14142320 -> 0.0054

XPL0

code Ran=1, CrLf=9;
code real RlOut=48;

func real MontePi(N);   \Calculate pi using Monte Carlo method
int  N;                 \number of randomly selected points
int  I, X, Y, C;
def  R = 10000;         \radius of circle
[C:= 0;                 \initialize count of points in circle
for I:= 0 to N-1 do
        [X:= Ran(R);
         Y:= Ran(R);
        if X*X + Y*Y <= R*R then C:= C+1;
        ];
return float(C)*4.0 / float(N);   \Acir/Asqr = pi*R^2/4*R^2 = pi/4
];

[RlOut(0, MontePi(        100));  CrLf(0);
 RlOut(0, MontePi(     10_000));  CrLf(0);
 RlOut(0, MontePi(  1_000_000));  CrLf(0);
 RlOut(0, MontePi(100_000_000));  CrLf(0);
]

Output:

    2.92000
    3.13200
    3.14375
    3.14192

zkl

fcn monty(n){
   inCircle:=0; 
   do(n){
      x:=(0.0).random(1); y:=(0.0).random(1);
      if(x*x + y*y < 1.0) inCircle+=1;
   }
   4.0*inCircle/n
}

Or, in a more functional style (using a reference for state info):

fcn monty(n){
   4.0 * (1).pump(n,Void,fcn(r){
      x:=(0.0).random(1); y:=(0.0).random(1);
      if(x*x + y*y < 1.0) r.inc(); 
      r
   }.fp(Ref(0)) ).value/n;
}

Output:

T(100,1000,10000,0d100_000,0d1_000_000,0d10_000_000)
   .apply2(fcn(n){"%10,d : %+f".fmt(n,monty(n)-(1.0).pi).println()})
       100 : -0.061593
     1,000 : +0.018407
    10,000 : -0.013993
   100,000 : -0.000833
 1,000,000 : -0.004385
10,000,000 : +0.000619